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ABSTRACT 


Integral  equations  (IE's)  are  widely  utilized  to  calculate  induced  currents  on  anten¬ 
nas  and  scatterers,  but  they  are  seriously  restricted  in  their  ability  to  handle  inhomoge¬ 
neous  penetrable  structures  having  multiwavelength  dimensions.  The  utilization  of  finite 
element  (FE)  techniques  has  not  been  as  pervasive  as  the  use  of  IE's.  The  IE  represen¬ 
tation  matrix  is  "full",  containing  few,  if  any,  zero  valued  elements.  The  techniques  for 
operating  on  these  large-sized  full  matrices  require  undesirable  amounts  of  processor 
time.  FE  techniques  produce  sparse  matrices  due  to  the  strictly  local  interactions  be¬ 
tween  discrete  unknowns.  The  application  of  FE's  to  unbounded  problems,  however, 
requires  supplemental^  enforcement  of  the  far-field  radiation  conditions.  The  Field 
Feedback  Formulation  (F^)  circumvents  the  full-matrix  computational  "bottleneck"  by- 
allowing  FE  based  numerical  methods  to  be  employed.  Even  though  the  resultant  sparse 
matrices  may  be  larger  than  the  "full"  matrices  discussed  earlier,  most  elements  have  a 
value  of  zero.  Numerical  procedures  exist  to  optimize  operations  with  these  sparse  ma¬ 
trices.  Calculational  speeds  can  be  orders  of  magnitude  faster.  Computer  techniques  to 
implement  and  validate  this  new  technique  are  the  basis  for  this  thesis.  Excellent 
agreement  with  classical  results  are  demonstrated. 
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I.  INTRODUCTION 


A.  HISTORY 

The  application  of  finite  element  techniques  to  evaluate  the  solution  of  differential 
equations  is  well  documented.  The  utilization  of  these  techniques  by  the 
electromagnetics  community  has  not  been  as  pervasive  as  the  use  of  integral  equations 
(IE's).  IE's  are  widely  utilized  to  calculate  induced  currents  on  antennas  and  scatterers, 
but  they  are  seriously  restricted  in  their  ability  to  handle  inhomogeneous  penetrable 
structures  having  multiwavelength  dimensions.  As  the  complexity  and  number  of  nodal 
degrees  of  freedom  grow,  the  size  and  dimension  of  the  representative  matrix  must  also 
grow.  This  matrix  is  "full ",  containing  few,  if  any,  zero-valued  elements.  Ihe  available 
numerical  techniques  for  operating  on  these  large-sized  full  matrices  require  undesirable 
amounts  of  processor  time. 

Differential  equation  (DE)  based  techniques,  such  as  the  finite  element  method, 
produce  sparse  matrices  due  to  the  strictly  local  interactions  between  discrete  unknowns 
which  result.  The  application  of  DE's  to  unbounded  problems,  such  as  those  of  scat¬ 
tering  and  radiation,  require  some  form  of  supplementary’  enforcement  of  the  proper 
far-field  conditions.  These  radiation  boundary  conditions  are  innately  incorporated  into 
integral  equations. 

B.  FIELD  FEEDBACK  FORMULATION 

The  Field  Feedback  Formulation  (F^)  circumvents  the  full-matrix  "bottleneck"  in  the 
computational  process  by  allowing  DE  based  numerical  methods  to  be  employed.  Even 
though  the  resultant  sparse  matrices  may  be  larger  than  the  "full"  matrices  discussed 
earlier,  most  elements  have  a  value  of  zero.  Numerical  procedures  exist  to  optimize 
operations  with  these  sparse  matrices.  Calculational  speeds  can  be  orders  of  magnitude 
faster  than  with  full  matrices  (Ref.  1].  Although  the  F* employs  sparse  matrices  to  rep¬ 
resent  the  fields  in  the  materials  being  considered,  it  does  require  augmentation  to  en¬ 
force  the  radiation  condition  at  infinity  on  the  scattered  fields.  This  comes  in  the  form 
of  a  feedback  matrix  composed  of  surface  integration  generated  elements.  A  concept 
evaluation,  for  a  special  axisymmetric  case,  was  already  accomplished,  as  detailed  in 
[Ref  2]  and  [Ref  3).  Computer  techniques  to  implement  and  validate  this  new-  technique 
are  the  basis  for  this  thesis. 
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C.  POTENTIAL  BENEFITS 

This  thesis  will  lead  to  an  increased  understanding  of  the  advantages  and  disadvan¬ 
tages  of  this  novel  computational  procedure  for  handling  geometrically  complex  material 
scatterers.  This  method  may  ultimately  allow  computer-aided  design  of  important 
electroma;;  .etic  structures  such  as  low-observable  aircraft,  high  efficiency  dielectric  lens 
antennas,  and  other  electromagnetic  scattering  occurrences  due  to  atmospheric  anoma¬ 
lies.  Structural  details  and  material  inhomogeneities,  as  well  as  physical  dimensions  (in 
multiple  wavelengths),  can  be  accommodated  using  the  Field  Feedback  Formulation. 
These  capabilities  far  surpass  those  which  are  possible  with  contemporarx-  integral 
equation  techniques  for  the  case  of  inhomogeneous  penetrable  scatterers  and  antennas. 
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II.  FORMULATION 


A.  INITIAL  NOMENCLATURE 

Assume  there  is  a  three  dimensional  object  that  is  infinite  in  one  direction.  Such  an 
object,  in  cross  section  could  look  like  Figure  I.  This  object  only  varies  in  two  dimen¬ 
sions,  and  therefore,  is  actually  a  two  dimensional  (2-D)  object. 


Figure  1. 


y 


The  wavenumber  Ag  is  defined  as 


,  2ff  2nfo 

where  is  the  free  space  wavelength  associated  with  an  electromagnetic  wave  of  fre¬ 
quency  and  c  is  the  speed  of  light.  The  X  and  Y  Cartesian  coordinates  are 
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wavenumber  normalized,  such  that  and  Y  =  k^.  This  coordinate  normalization 

will  be  used  throughout  this  development.  A  similar  technique  could  also  be  developed 
in  the  polar  coordinate  system.  The  magnetic  field  will  also  be  normalized  such  that 
tt  =  .  where  »/o  =  =  I20n  »  377^2,  is  the  impedance  of  free  space  and 

is  the  usual  magnetic  field  in  units  of  A/m.  Thus  the  normalized  H  has  the  same  V/m 
units  as  E.  Potentials  may  be  defined  as, 

E,ix^)=UX,Y)  (2.1) 

for  the  transverse  magnetic  (TM)  case,  with  E„  E,,  and  //,  =  0  ,  and 

(2.2) 

for  the  transverse  electric  (TE)  case,  with  //„  //,  and  £,  =  0  . 

Our  objective  is  to  calculate  the  scattered  fields  for  an  arbitrary  (2-D)  penetrable 
object  using  the  Field  Feedback  Formulation  (F  ).  As  shown  in  Figure  2,  a  familiar 
closed  loop  system  illustrates  the  relationship  between  the  incident,  scattered,  total  and 
far  fields. 


Figure  2.  Field  Feedback  Formulation 

At  point  1,  the  incident  field  drives  the  total  system.  The  incident  field  at  1,  combined 
with  the  scattered  field  at  4,  forms  the  total  field  at  2.  This  total  field  forms  the 
boundary  conditions  that  drive  the  finite  element  program  at  2.  At  point  2,  initially,  the 
incident  field  drives  the  U  operator.  U  represents  the  feed  forward  operator  in  the  f! 
This  operator  uses  a  finite  element  technique  to  solve  the  boundary  value  problem.  At 
"oint  3,  the  boundary  conditions  are  solved  for  the  object  perimeter  potentials  and  the 
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normal  derivative  of  these  potentials.  T  represents  the  field  feedback  operator  that  takes 
the  perimeter  potentials  and  associated  derivatives  and  provides  the  scattered  fields  at 
point  4,  the  offset  boundar}'.  The  fields  at  point  1  and  4  are  added  (unlike  the  familiar 
feedback  or  control  system  where  negative  feedback  is  employed).  At  point  2,  a  com¬ 
bined  incident  and  scattered  field  exists.  These  fields  are  the  combined  boundary  con¬ 
ditions  for  the  finite  element  boundary  value  program.  These  combined  fields  (total 
fields)  are  then  applied  to  the  U  operator  to  calculate  the  perimeter  fields  and  the  asso¬ 
ciated  derivatives  on  the  objects  perimeter.  This  looping  may  be  repeated  until  a  steady 
state  condition,  at  point  3,  is  reached.  The  existence  of  a  steady  state  condition  assumes 
stability.  Stability  for  physical  systems  should  not  be  a  problem.  However,  when 
mathematically  modeled,  instabilities  may  result.  The  error  magnification  or  condition 
number  of  the  system  must  also  be  seriously  considered  if  this  iterative  looping  process 
is  to  be  used.  The  alternative  approach  is  to  form  an  equivalent  system,  where: 

equivalent  operator  =  U •[_!  —  T •  TT"', 


with 


I  —  Identity  Matrix 
and 

T*  U  =  Combined  Effects  of  the  T  and  U  operators . 

Either  approach  is  viable  since, 

^  total  =  ^Incident  +  T  •  V  •  ij/ incident  +  [T  •  Vf  *  ^Incident  +  •  •  •  =  [/  —  T»  T’]  incident 

This  becomes  more  obvious  when  is  factored  from  the  equality  leaving, 

I  4-  r.  U+(T^bf  +  ,  .  . 

Placing  this  in  closed  form. 


^(r.  C/f  = 

«=o 


1 

I-T^U 


=  il- 


r.  CO"’. 
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Finding  the  equivalent  operator  will  require  a  matrix  inversion  and  for  very  large  prob¬ 
lems  this  may  lead  to  excessive  computation  times.  The  matrix  inversion  technique  will 
be  investigated.  The  fields  may  then  be  extended  to  any  point  in  space  using  a  far  field 
Green  s  function  surface  contour  integral.  This  far  field  pattern  is  available  at  point  5. 

The  incident  field  is  usually  produced  by  a  plane  wave  generator.  This  provides  the 
boundary  conditions  on  the  offset  contour.  This  contour  is  called  "offset"  since  it  is 
approximately  the  same  "shape"  as  the  objects  perimeter  but  is  slightly  larger.  The  dis¬ 
tance  between  the  perimeter  and  this  contour  is  called  the  offset  distance  and  will  be 
discussed  in  Chapter  III.  The  boundaty  conditions  may  be  any  desired  field  or  wave  that 
satisfies  Maxwell's  equations.  These  waves  may  arrive  from  any  direction  and  be  of  any 
magnitude.  The  boundary  conditions  may  also  be  a  composite  of  any  number  of  waves 
since  superposition  does  apply  to  these  systems.  A  user  provided  subroutine  is  necessaiy 
if  conditions  other  than  a  single  plane  wave,  cylindrical  mode  (of  arbitraiy  mode  num¬ 
ber)  or  individual  input  boundaiy  condition  is  desired. 

B.  MAXWELL'S  EQUATIONS 

Maxwell's  equations  can  be  written  using  our  previous  normalizations  as, 

VxE^fiJI  (2.3) 

and 

V  X  /7  =  c,E.  (2.4) 

d 

[Ref.  4  ]  Let  ,  with  similar  definitions  for  Dy  and  Dy.  Equations  2.3  and  2.4 

can  be  further  expanded  into  the  0^,  Dy  and  Dy  components  such  that. 


(2.5) 

HyHy  =  -DxEy 

(2.6) 

(IfHy  =  DxEy  —  DyEx 

(2.7) 

CyEy.  =  DyHy 

(2.8) 

lyEy  =  ~DxHy 

(2.9) 

tfEy  =■  DxHy  —DyHje- 

(2.10) 
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For  the  TM  case,  with  propagation  in  the  z  direction, 


DyE, 

JJ  = — L_L 
"jc  Mr 


and 


-DyE, 


Note  that  the  H,  field  =  0.  These  two  equations  can  be  combined  to  form, 


X  2". 


(2.11) 


Similarly  for  the  TE  case,  with  propagation  in  the  z  direction, 


£  =  Xv^2  X  2. 


(2.12) 


Substituting  equations  2.5  and  2.6  into  equation  2.10  yields, 

■DxE, 


^rE,  =  Dx\ 


Mr 


(2.13) 


Substituting  equation  2.1  into  equation  2.13  yields, 


Er'Ai  +  l/lij  +  Dy^  (/r, 


=  0. 


(2.14) 


Equation  2.14  can  be  further  simplified  to, 

V.[-;j;:V^,]+a,=0.  (2.15) 

Similarly,  equations  2.2,  2.8,  2.9  may  be  substituted  into  equation  2.7.  This  yields, 

V.[-^V|/r2]+Mr«/'2*0.  (2.16) 

Equations  2.15  and  2.16  are  T.M  and  TE  duals.  These  two  differential  equations  describe 
the  potentials  inside  the  object  of  interest.  Defining  “  =  «  and  c,  =  fi  for  the  TM  case 
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and  =  a  and  fi=P  for  the  TE  case  and  substituting  these  new  definitions  into 
equations  2.15  and  2.16  yields  one  differential  equation, 

V .  [aV^fr]  + =  0.  (2.17) 

C.  VA.RUTIONAL  EQUIVALENCE  TO  THE  DIFFERENTIAL  EQUATION 

The  Euler-Lagrange  variational  formulation  is  based  on  the  stationarity  of  a  func¬ 
tional,  of  the  function  \jj  and  its  first  derivatives.  (Ref  4  ) 


/  =  JJ  FiX,Y,^,Vil/)dXdY  (2.18) 

inside  S 

It  can  be  shown  that  the  first  variation  of  the  functional  is  zero,  SI=0  ,  if  the 
Lagrangian,  F,  satisfies  the  Euler  -  Lagrange  equation: 


e  (  dF  \  c  (  cF 

cX  \  d{Dx4>)  J  5K  V  SiDyxlj) 


(2.19) 


The  problem  thus  becomes  to  find  the  F,  which  when  substituted  into  equation  2.19. 
yields  the  original  differential  equation.  The  Lagrangian, 


can  be  simplified  to. 


F=aV.A.V^-/?^l 


When  equation  2.20  is  substituted  into  equation  2.19, 


dF 


d{Dx^) 


cF 

diDy^k) 


=  2a  Dy^ 


(2.20) 


Therefore, 


{2<xDx4>)  +  -fy  OaDyip)  +  2)5^  =  0 
or. 


which  simplifies  to, 


=  0 


V  •  [aVi/f]  +  /?^  =  0. 


(2.21) 


Therefore,  the  general  functional  has  been  found  since  equation  2.21  and  equation  2.17 
are  identical.  Substituting  the  a  and  /?  definitions  into  equation  2.20  yields, 


(2.22) 

and 

(2.23) 


Equations  2.22  and  2.23  are  integrated  over  the  interior  to  S  with  known  boundarv- 
conditions  for  either  or  on  S.  To  physically  interpret  the  variational  formulation, 
it  is  noted  that, 


and 

Therefore, 


t,E  •  E  dXdY 


s 


and 


9 


-JJ. 


i2-J  ^  •  E  ~  tirH  •  H  dX  dY. 

s 


Substituting  Maxwell's  equations  into  /,  gives. 


£).//-(V  X  H)‘EdXdY 


"■f! 


/,  =  V.(£  X  H)dXdY 


/,  =  E  X  H»ndL 

^dS 


Note  that  E  x  H  is  the  complex  oscillator}'  Poynting  vector.  It  is  different  from  the 
usual  Poynting  vector  of  £  x  H\  Thus,  both  functionals,  /,  and  are  proportional  to 
the  complex  phasors  for  oscillator}'  power.  The  oscillatory  power  is  the  phasor  repres¬ 
enting  the  excess  of  instantaneous  radiated  power  minus  the  average  radiated  power. 

D.  FINITE  ELEMENT  BOUNDARY  VALUE  SOLUTION 

With  the  discussion  of  the  variational  mathematics  completed,  a  simple  rectangular 
boundary'  value  problem  will  be  discussed  in  detail.  The  rectangular  geometry'  allows  for 
an  easier  formulation  but  in  no  way  limits  the  solution  from  being  extended  to  more 
complicated  object  geometries.  Figure  3  on  page  1 1  shows  the  region  of  concern. 
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y  t 


0 

0  Q  X 

Figure  3.  Rectangular  Object 

Consider  spanning  the  rectangular  region  by  a  triangular  mesh,  as  shown  in 
Figure  4  on  page  12.  Both  a  and  may  be  functions  of  position  but  may  not  vary 
within  an  individual  triangular  element.  Given  a  fine  enough  mesh  structure,  a  smooth 
transition  in  material  properties  may  be  approximated.  For  notational  simplicity  this 
positional  dcpendance  will  not  be  carried  forward.  The  variational  approach  will  yield 
the  solution  to  the  boundarj'  value  problem  by  finding  the  \l/{x,y)  which  gives  the  sta¬ 
tionary  value  of, 

nb 

(aV(^  •  dx  dy 

which  is  constrained  on  the  boundary  by  the  previously  specified  boundary  conditions. 
(Ref.  5  J 
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r 


1-01  M  M+1 

Figure  4.  Rectangular  Mesh 


1  he  values  of  rif  at  the  interior  nodes  become  the  discretized  unknowns: 

=  ^{Xi,  Yj)  I  —  1. . .  M  andj  -  I . ..  N 


where 

X,  =  i^AX 


and 


for  the  uniform  mesh  structure  W'ith  AX 


(Af+  1) 

W+J  A-+I 


and  AY  = 


(/V+1) 


.  Approximating, 


^ix,y)=  Yj 

M  JmO 


which  includes  the  known  boundary  nodal  values  and  ^u-nj  for  j  =  0  and  N  +  1, 
and  a  and  for  i  =  0  and  M  +  1,  linear  pyramidal  basis  functions,  j;),  which 
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have  unit  value  at  the  (i,j)  node  and  zero  value  at  all  surrounding  nodes.  Figure  5(a)  is 
a  top  view  of  the  pyramidal  basis  function,  while  Figure  5(b)  is  a  perspective  view. 


Figure  5.  Pyramidal  Basis  Function,  (a)  top  view,  (b)  perspective  view 


The  functional  I  w'ill  thus  be  a  discrete  function  of  each  of  the  nodal  values  of  r{/, 

i  ^0,  i»  •••*  ^+l)• 

The  approximate  discrete  solution  will  be  found  by  the  system, 

~ —  =  0,  for  m  =  I . . .  Af  and  n  =  I  ...  N. 


Now, 


dl 


na 


rb 


*'0  •'o 


oVt/f 

^^m,n 


S^m.n 


dxdy^a 
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where, 


M+\  N+\ 


1=0  y=o 


The  gradient  of  the  basis  function  is, 


^«m,  «(•*■.>')-  Q. 

^Tm,n 


and  the  basis  function  is. 


^rm,  n 

Therefore,  the  system  of  equations  to  solve  becomes. 


nb 

(aVu^  „  •  -  /?«„,  „^)  dr  ^y  =  0,  for  /«  =  1  . . .  M  and  «  =  1  . . .  A’. 

) 


Substituting  for  ViJ/  and «//  in  terms  of  the  nodal  values  of  \l/  for 
m=  \ . . .  M  and  «  =  1  . . .  A"  gives. 


W+'  v+i 

V— I  r“  r* 


i=l  y=0 


0  •'0 


where, 


ra  r* 

(aV«„.  «•  t/jf  </>•  =  f{(m, «),(/, y)}. 

0  •’0 


Regrouping  this  to  put  the  known  nodal  values  on  the  right  hand  side  gives  for 
m=  \  ...  M and n  =  1 . . .  A' (interior  nodes). 


M  .V 

X  «).(', 7’)]  =  -  'Yj 

1=1  y=l  Boundary 

Nodes  Only 
for(ij) 
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By  renumbering  the  nodes  using  a  single  index  for  the  unknown  nodal  values, 
k  —  N{i  —  l)+J  and  I  =  Xim  —  1)  +  «  and, 

;t=i  Jt 

The  functional  F{1,  k)  =  0  if  nodes  /  and  k  are  not  both  associated  with  at  least  one 
common  triangular  element.  Therefore,  and  U/U/,  will  be  zero  except  in  triangles 

where  /  and  k  both  appear  as  nodes.  For  the  example  mesh  structure  of  Figure  3  on 
page  II,  this  produces  a  banded  matrix,  F.  This  sparse  matrix  can  be  easily  displayed 
by  first  defining  column  vectors  of  the  unknown  nodal  values  in  the  mesh 

%  =  C'A/,  li 2«  •••»  •A/,  v]  • 

The  F  -  matrix  elements  F  [(m,  «),(/,y)3ire  zero  unless  the  (m.n)  node  shares  at  least  one 
element  with  the  (i,j)  node.  Thus,  the  node  values  in  \p,  will  be  coupled  only  to 
lA/.i,  iA,_i  and  any  associated  boundary  nodes.  This  is  written  as, 

[/i JV,  +  [B,]?,  + 

where  A„  B,  and  are  A'  x  A’  complex  arrays  and  P,  is  a  A'  x  A’j,^  array  where  A'i,_  is  the 
number  of  boundary  nodes  associated  with  the  i-th  column.  This  appears  as. 


If  we  denote  as  the  initial  boundary  values  and  as  the  final  boundary  values 
then  becomes  only  the  boundaiy  conditions  on  the  top  and  bottom  at  y  #  0  and 
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iV  +  1.  Note  that  the  above  system  is  tri-block  in  nature  and  has  a  large  number  of  zero 
valued  elements.  The  zero  valued  elements  were  omitted  for  clarity.  Each  element  is 
actually  a  matrix  and  therefore  it  is  evident  how  sparse  this  system  is.  Each  of  tlie 
Bi,  C,  and  P,  matrices  are  equally  sparse,  however,  a  global  symmetry  does  not  ap¬ 
pear. 

E.  EVALUATION  OF  THE  F  -  MATRIX  CONTRIBUTIONS 

Given  an  arbitrary  element  as  shown  in  Figure  6,  the  potential,  i/r,  can  be  linearly 
approximated  by  [Ref.  5  ], 

.'^3. 

3 

k=\ 

1 


3 


2 

Figure  6.  Mesh  Element 

where  is  the  nodal  value  at  the  k-th  node,  [T]  *3x3  Transform  Array 

and 
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Ui/^x,y)  =  Linear  basis  functions  for  the  k-th  node 


=  {x,y,  1) 


T'2.k 

h.u 


-  ^l,k^  +  ^2,ky+  T^3,k- 


Note  that, 


1  ,  fork  =  nx\ 
0  ,  for  k  i=-  mj 


It  can  be  shown  that, 


CV2-^3) 

(•>(■3  -  -^2) 

(•’!'2>'3  ~  •’^iV'2) 


0  '3  ~>  l) 

(Jf]  -  -^3) 

(jfiV,  -  Jft^-3) 


iy\  ">'2) 
U2  - 

(jfl>2  -  X2.V^) 


where  M 1  =  triangle  area,  and 


2A  =  det 


y\ 

Xl  >2 
'’^3  >^3 


1 

1 

1 


2 A  =  (;CiV3  +  jr3^v,  +  x^y^)  -  +  x^y^  +  jr^y,) 


Furthermore,  the  linear  basis  functions,  can  be  interpreted  as  relative  areas  of 

the  triangle  shown  in  Figure  7  on  page  18. 
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1 


Figure  7.  Basis  Functions  Interpretation 


A  =  A^  +  A2  +  A^ 


is  constant  as  varies  and, 


Ukix,y) 


A 


Within  a  given  triangular  element,  the  evaluation  for  k  =  1,  2,  3  and  /  =  1,  2,  3  in  the 
element  is  of  interest  and, 


Vu,  •  Vufc  -  dx  dy. 


inside 


triangle 


q 

The  assumption  is  made  that  a  and  P  will  be  approximated  as  constants  within  the  tri¬ 
angle.  These  material  constants  can,  however,  vary  from  element  to  element.  Taking 
the  gradient  terms  first, 


Vu/«Vujt==(r,j^»  Ti  ,A-  Tj  /). 
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Therefore, 


j  j 

triangle 

==  ^2.  *^2,/)  I -^(7 1* 

where  A,  is  the  area  of  the  q”'  element.  Next  it  can  be  derived  that, 


r-jy  Ml  y  fork  ¥=h 
uiu^dxdy  =  ^ 

I  -^Ml  ,  forked 

triangle 

More  generally,  with  each  u,  raised  to  an  integer  pow'er,  n„ 


f  f u”'u2^U3^dxdy  =  2 M 1 


triangle 


(rt,  +  /J2  +  ^3  +  2)! 


The  final  result  is. 


/:)  =  JJ  (aVu,.  -  ^uiuf}dxdy 

triangle 

q 

F.  GREEN'S  FUNCTION  CONTOUR  INTEGRAL 

The  scattered  fields,  i/t,  from  an  arbitrary  object  in  a  vacuum  satisfying  Helmholtz's 
equation  (see  Figure  8  on  page  20  ), 

V^i/i  +  AV  =  0, 

are  (Ref  6  ], 
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contour 


Figure  8.  Green's  Function  Integration 


where  the  Green's  function  is, 


^'1'  A  ^ 


=  n*V^  on  the  contour 


-Ci.i  ihl.  t,  1  r  r-  h 
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The  Hankel  functions  of  the  second  kind  of  order  zero  and  one,  and  ,  will 

present  a  problem  for  the  numeric  integration  discussed  in  Chapter  V.  The  imaginary 

portion  of  these  functions  rapidly  approaches  negative  infinity  as  the  argument  ap- 

dif/ 

proaches  zero.  The  is  obtained  by  a  finite  difference  method  using  the  field 
boundary  conditions  on  surface  B  (boundary  conditions)  and  the  calculated  field  condi¬ 
tions  on  surface  P  (object  perimeter).  This  results  in, 

Oij/  ^boundary  ^perimeter 

dn  offset  distance 


It  •will  be  sho-wn  that  to  maximize  the  accuracy  of  the  Green's  function  the  offset 

dij/ 

distance  should  be  made  as  large  as  possible.  This,  however,  causes  the  to  be  in¬ 
ert 

accurate.  Thus,  an  optimal  condition  must  be  found  that  maximizes  the  accuracy  of  the 


entire  numeric  integration.  Such  a  condition  does  not  maximize  the  accuracy  of  any  one 


of  the  contributing  parts  to  the  Green's  function  integrand. 


G.  FAR-FIELD  EVALUATION 

When  the  Green's  function  integral  discussed  above  is  used  for  far  field  calculations 
several  simplifying  relationships  develop.  These  simplifications  require  a  less  demanding 
numerical  integration.  To  be  in  the  far  field  region  three  conditions  must  exist. 


I F I  >  Z),  I F I  >  /(,  and 


F  > 


2D^ 


where  /(,  is  the  free  space  wavelength  and  D  is  the  maximum  dimension  of  the  object. 
[Ref  7  ]  As  X  approaches  infinity. 


and 


This  requires. 


G(F|F') 


and 
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cG  ^  D  ^0 
=  -  n  •  R—r 
cn  4 


In  the  far  field 


2j 


nkc.R 


,-yM 


yj and  n  •  R  -*  n  •  r.  Thus, 

g-JkoR  ,  gJI^cT  • 


Therefore, 


G(r  Ip)  = 


J_ 

4 


and 


With  these  new  definitions  substituted  into  the  original  Green's  function  integral 
equation. 


^ scaitereJ<^)  ~ 


^nk^r 


,-yv 


J-^  +  kQn*r\li{r  ) 
C  fl 


*'C 

Note  that  this  equation  is  partitioned  into  a  distance  dependent  term  and  a  theta  depend 
term.  The  theta  dependent  term  may  be  defined  as 


/  = 


Bn 


+  A:o«-r^(P)p^ 


cos  0 


e/c'. 


•'C 

The  two  dimensional  bistatic  radar  cross  section  (RCS)  per  unit  length  of  the  cylindrical 
structure  may  now  be  defined  as. 


RCS  =  o(0",  </>')  =  lim 

'  '  r-»oo 


27rrR^ 


=  lim 

/"-♦OC 


27rr| 

W\^ 
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a  —  lim  litr  =  — rr 

'■-oo  %nhr  Ak(s 


where  the  wavenumber,  and  the  incident  field  is  assumed  to  be  of  unit 

^*0 

magnitude,  =  1.0. 
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III.  MESH  GENERATION 


A.  INTRODUCTION 

A  display  or  plot  of  the  computer  generated  mesh  structure  is  not  required  for  the 
problem  solution,  however,  it  provides  an  immediate  visual  confirmation  that  the  in¬ 
tended  problem  geometry  has  been  entered  correctly  into  the  computer.  A  large  amount 
of  initialization  data  is  required  for  even  the  most  rudimentary’  problem.  For  this  reason, 
the  input  data  is  provided  to  the  mesh  generation  program  via  a  data  file.  Modifications 
are  possible  at  a  later  time  wth  minimum  effort. 

B.  INPUT  DATA 

The  input  data  file  is  called  INPUT.DAT  and  contains  25  input  fields.  Of  these 
fields,  14  relate  directly  to  the  generation  or  display  of  the  mesh  structure.  Only  the 
object  surface  coordinates  require  adherence  to  a  specific  format.  All  other  data  need 
only  be  of  the  correct  type  (i.e.,  character,  integer  or  real).  A  brief  description  of  each 
field  is  provided  below. 

•  Field  I  is  a  label  of  no  more  than  12  characters.  This  label  is  for  the  plot  and  input 
data  file. 

•  Field  2  is  a  character  flag  that  specifies  the  input  coordinate  system.  If  set  to  "R" 
or  "r“,  the  rectangular  system  is  used.  If  set  to  "P"  or  "p",  the  polar  system  is  used. 
If  a  "P","p"  ,"R"  or  "r"  is  not  detected,  an  error  is  returned  to  the  display. 

•  Field  3  is  a  character  flag  that  if  set  to  “I"  or  "i”  will  cause  several  intermediate 
values  to  be  stored  to  disk  during  the  Finite  Element  Boundary  Value  (FEBV) 
program  execution.  This  option  was  used  during  the  debug  process. 

•  Field  4  is  a  character  flag  that  if  set  to  "D"  or  ~d"  will  create  a  DISSPL.A 
FORTRAN  program  capable  of  replicating  the  input  object,  in  wavenumber  nor¬ 
malized  coordinates,  on  mainframe  computers  having  a  DISSPLA  graphics  pack¬ 
age.  DISSPLA  is  a  subroutine-based  language.  The  generated  program,  called 
DISSPLA. FOR,  is  a  compilation  of  four  subroutines  calls  per  element.  This  file 
can  get  very'  large  for  dense  mesh  structures. 

•  Field  5  is  a  character  flag  that  if  set  to  "U”  or  "u"  will  cause  a  uniform  material  to 
be  assumed.  In  the  uniform  case,  no  material  interface  exists.  This  option  was 
only  used  to  verify  the  Finite  Element  solution  accuracy. 

•  Field  6  is  a  character  flag  that  if  set  to  ".M"  or  "m"  will  cause  only  the  mesh  to  be 
generated.  This  is  very  useful  when  first  starting  a  problem  and  the  optimal  mesh 
structure  has  not  been  determined. 

•  Field  7  is  a  real  number  specifying  the  desired  mesh  resolution  in  wavelengths.  The 
mesh  resolution  determines  the  dimension  of  the  mesh  elements. 


24 


•  Field  8  is  a  real  number  specifying  the  distance,  in  wavelengths,  between  the  object 
perimeter  and  the  offset  boundaiy  contour. 

•  Field  9  is  a  real  number  that  specifies  a  multiplicative  scaling  factor  for  the  nu¬ 
merical  integration  stepping  function  discussed  in  Chapter  V. 

•  Field  10  is  a  bias  term  that  can  be  used  to  shift  the  numerical  integration  stepping 
function.  This  term  is  used  for  distances  less  than  1.0. 

•  Field  1 1  is  a  bias  term  that  can  be  used  to  shift  the  numerical  integration  stepping 
function.  This  term  is  used  for  distances  greater  than  1 .0. 

•  Field  1 2  is  a  real  number  specifying  the  maximum  distance  beyond  which  no  further 
contribution  to  the  Green's  Function  Integral  is  made.  If  this  feature  is  not  de¬ 
sired,  this  term  should  be  made  larger  than  the  objects  maximum  dimension  plus 
twice  the  offset  distance. 

•  Field  13  is  an  integer  specifying  the  number  of  input  data  points. 

•  Field  14  is  an  integer  specifying  the  angular  resolution,  in  degrees,  desired  for  the 
final  radar  cross  section  calculation. 

•  Field  15  is  an  integer  specifying  the  mesh  generation  technique. 

•  Field  16  is  an  integer  specifying  the  perimeter  node  from  which  the  bisection  seg¬ 
ment  originates.  This  node  is  called  the  "start  node". 

•  Field  17  is  an  integer  specifying  the  perimeter  node  on  which  the  bisection  segment 
terminates.  This  node  is  called  the  "stop  node". 

•  Field  18  is  a  pair  of  real  numbers  (on  two  lines)  specifying  the  x  and  y  coordinates 
by  which  the  object  will  be  displaced. 

•  Field  19  is  a  real  number  specifying,  in  wavelengths,  the  desired  distance  between 
the  origin  and  the  first  input  data  point.  Fields  18  and  19  when  used  together,  al¬ 
low  an  object  to  be  placed  at  any  position  and  scaled  to  any  size. 

•  Field  20  is  a  pair  of  real  numbers  (on  tv/o  lines)  specifying  the  real  and  imaginaiy- 
parts  of  ~  for  the  TM  case  or  —  for  the  TE  case, 

•  Field  21  is  a  pair  of  real  numbers  (on  two  lines)  specifying  the  real  and  imaginarj- 
parts  of  e,  for  the  T.M  case  or  n,  for  the  TE  case. 

•  Field  22  is  a  character  flag  that  if  set  to  "P"  or  "p"  enables  a  plane  wave  generator. 
The  plane  wave  is  propagating  down  the  y  axis,  and  generates  an  condition 
on  the  offset  boundary  contour.  If  the  flag  is  set  to  "C"  or  "c",  a  cylindrical  mode 
generator  is  enabled.  This  generates  an  Ex)  cos  n(t>  condition  on  the  offset  boundary 
contour.  If  a  "P"  ,  "p",  "C"  or  "c"  is  not  detected,  then  manually  input  boundary 
conditions  must  follow,  and  fields  23  and  24  are  not  used. 

•  Field  23  is  a  real  number  specifying  the  wave  amplitude. 

•  Field  24a  is  a  real  number  specifying  the  wave  frequency  (in  Hz). 

•  Field  24b  (only  for  cylindrical  case)  is  an  integer  specifying  the  mode  number,  n. 

•  Field  25  is  the  object  perimeter  data,  in  either  polar  or  rectangular  form. 
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An  example  of  an  input  data  file  for  a  homogeneous  circular  cylinder  is  provided  in 
Appendix  A.  The  majority  of  the  input  data  is  echoed  to  the  computer  display  and  a 
system  "pause"  is  initiated  to  allow  for  user  inspection.  The  program  may  be  aborted 
or  continued  at  this  time.  The  initial  object  dimensionalization  has  already  occurred, 
and  the  number  of  unknowns  and  the  maximum  unknown  width  is  displayed.  These 
factors  give  an  excellent  indication  of  the  expected  run  time  for  the  FEBV  routines.  For 
example,  a  problem  with  512  unknowns  and  a  maximum  unknown  width  of  31  took  806 
seconds  to  execute  while  a  run  w-ith  8  unknowns  and  a  maximum  unknown  width  of  3 
took  only  10  seconds  to  execute.  These  times  are  for  a  Intel  80386  based  personal 
computer  with  an  Intel  80287  co-processor  chip  calculating  the  fields  for  a  circular  cyl¬ 
inder.  The  FEBV  routines  are  the  ne.xt  code  block  to  execute  after  the  pause  is  cleared. 

C.  MESH  GENERATION  PROGRAM 

The  mesh  generation  program  consists  of  seven  subroutines.  These  subroutines  are 
an  integral  part  of  the  finite  element  program  and,  therefore,  were  not  separated.  These 
routines  are  discussed  below'. 

1.  lO  (Input/Output) 

This  subroutine  reads  the  information  contained  in  the  INPUT.DAT  file  dis¬ 
cussed  earlier.  A  two  dimensional  object  can  be  described  in  any  number  of  w'ays, 
however,  for  simplicity,  the  polar  and  rectangular  coordinate  systems  are  used.  In  either 
case,  the  initial  assumption  is  that  all  data  points  are  referenced  to  a  local  origin.  This 
local  origin  can  be  offset  by  any  desired  amount  using  field  18.  This  offset  is  independ¬ 
ent  of  the  entered  data  points  or  any  size  scaling  provided  by  field  19.  A  plot  label  file, 
named  TEXT.LBL,  is  created  and  the  initial  object  perimeter  (coded  for  a  display  in 
blue)  is  written  to  the  output  file,  PLT.DAT.  These  data  points  describe  the  perimeter 
of  the  object.  This  will  later  prove  helpful  in  determining  the  conformity  of  the  gener¬ 
ated  mesh  to  the  input  perimeter.  All  subsequent  screen  w'rites  are  coded  for  a  display 
in  green.  Two  example  objects  are  showm  in  Figure  9  on  page  27.  These  objects  will 
be  used  throughout  this  chapter.  The  circular  cylinder  was  generated  by  a  separate 
computer  program.  The  "horseshoe"  shaped  object  w'as  manually  input.  Graph  paper 
was  used  to  determine  the  x  and  y  coordinates  of  the  28  unequally  spaced  perimeter 
points. 


2.  Rotate 

Tliis  subroutine  reorders  the  input  data  points  to  allow  for  any  desired  bisection 
segment  start  and/or  stop  node.  This  subroutine  is  only  used  if  manual  selection  of  the 
bisection  segment  start  and/or  stop  nodes  is  requested. 

3.  Bound  (Boundary) 

This  subroutine  sub-divides  the  object  perimeter  based  on  the  mesh  resolution 
specified  in  Field  7.  The  mesh  resolution  is  the  approximate  length,  specified  in  wave¬ 
lengths,  that  the  user  desires  the  perimeter  to  be  divided  into.  The  division  of  the  per¬ 
imeter  is  based  on  linear  interpolation  between  input  data  points.  A  new  perimeter  node 
is  placed  at  each  of  the  sub-division  points.  Additional  input  data  points  are  necessary- 
in  areas  of  rapid  change  to  allow  for  a  correct  object  perimeter  representation.  1  he  ob¬ 
ject  bisection  e.xtends  from  the  "start  node’  to  the  "stop  node".  These  nodes  arc  user 
specified  in  Fields  16  and  17  and,  in  general,  divide  the  object  in  half.  The  bisection 
should  be  arranged  such  that  the  object  width,  perpendicular  to  the  bisection  segment, 
is  minimized.  Thus,  a  long  slender  object  should  be  oriented  for  a  major  axis  bisection. 
Whether  the  major  axis  is  oriented  vertically  or  horizontally  is  of  no  importance.  This 
is  demonstrated  in  Figure  10  on  page  28  (right  side). 
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BISECTION  SEGMENT 


Figure  10. 


“OBJECT  PERIMETER 

Bisection  Segment  and  Row  Arrangement 


For  more  complex  objects,  the  bisection  is  not  a  straight  line,  but  rather  a  scries 
of  line  segments  that  approximate  a  curve.  The  number  of  perimeter  nodes  on  the  left 
side  of  the  bisection  segment  must  equal  the  number  of  perimeter  nodes  on  the  right  side 
of  the  bisection  segment.  The  bisection  segment  also  contains  this  same  number  of 
nodes.  Thus,  the  program  must  adjust  the  requested  mesh  resolution  to  ensure  proper 
nodal  spacing.  This  allows  for  a  piecewise  continuous  segment  to  cross  the  object. 

4.  Normal 

This  subroutine  calculates  the  unit  normals  to  the  object  perimeter  at  each 
newly  established  perimeter  node.  The  normal  is  a  perpendicular  constructed  to  the 
chord  connecting  the  two  nodes  adjacent  to  the  node  for  which  the  calculation  is  oc¬ 
curring.  This  perpendicular  originates  at  the  current  node  and  is  of  unit  length.  See 
Figure  1 1  on  page  29. 
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5.  Nodset  (Node  Set) 

This  subroutine  uses  a  two  sweep  technique  to  compute  the  number  of  nodes 
on  each  nodal  row.  The  rows  are  labeled  I  =  1,  2,  3,  1  maximum  (IMX)  consec¬ 

utively  from  the  top  of  the  object  to  the  bottom.  Segments  normal  to  the  perimeter 
surface,  calculated  in  NORMAL,  connect  the  perimeter  nodes  to  the  offset  boundary 
contour.  Two  such  completed  normal  segments,  one  on  each  end,  complete  a  nodal  row. 
When  all  rows  are  complete,  a  second  contour  has  been  established  that  approximates 
the  object's  perimeter.  This  contour  is  displaced  from  the  perimeter  by  the  distance 
specified  in  Field  8.  This  contour  provides  the  boundary  conditions  for  the  FEBV  rou¬ 
tines.  The  optimal  selection  of  this  offset  distance  will  be  a  topic  of  Chapter  IV. 

With  the  object  divided  into  rows,  each  row  can  be  further  divided  into  equally 
spaced  nodes.  Based  on  the  requested  resolution  and  whether  the  objects  dimensions 
are  expanding  or  contracting,  the  nodal  spacing  is  adjusted  to  keep  the  elements 
approximately  the  same  size.  Two  adjacent  nodal  rows  form  an  elemental  row.  These 
rows  are  given  the  same  label,  I  =  1,  2,  3, ...,  IMX  as  the  upper  nodal  row.  The  nodes 
on  an  elemental  row  are  numbered  consecutively,  starting  with  the  left-most  node.  1  his 
node  is  always  an  offset  boundary  contour  node.  When  the  upper  nodal  row  is 


numbered,  the  process  is  continued  for  the  lower  nodal  row.  An  example  of  this  element 
row  numbering  scheme  is  shown  in  Figure  1 2  on  page  30. 


Figure  12.  Element  Row  Numbering  Scheme 


A  mesh  orientation  attribute  is  set  for  the  left  and  right  portion  of  each  element  row. 
This  attribute  determines  if  the  object  has  a  clockwise  or  counter-clockwise  orientation. 
See  Figure  13  on  page  31.  Switching  between  the  four  available  mesh  orientations  al¬ 
lows  for  a  mesh  that  more  accurately  conforms  to  the  input  object  perimeter  without 
resulting  in  a  disproportionate  mesh  structure. 


There  is  an  additional  requirement  that  the  number  of  nodes  on  a  given  left  or  right  half 
row  must  be  only  one  more,  equal  to,  or  one  less  than  the  number  of  nodes  on  the  ad¬ 
jacent  left  or  right  half  row.  Thus,  a  two  sweep  process  is  utilized.  This  process  ensures 
that  this  requirement  is  met  and  that  the  first  and  last  row  have  only  two  nodes.  The 
second  and  second  from  last  row  (if  present)  must  have  three  nodes.  It  is  possible,  as 
shown  in  Figure  14  on  page  32,  to  generate  a  mesh  that  has  only  three  rows.  This  object 
was  input  as  a  circular  cylinder.  It  is  evident  that,  due  to  a  small  number  of  elements, 
the  generated  mesh  more  accurately  resembles  a  square.  This  will  be  a  problem  for  cal¬ 
culating  the  scattered  fields,  however,  the  internal  fields  can  be  accurately  approximated. 
In  this  special  case,  the  second  and  second  from  last  rows  are  the  same.  It  will  be  shown, 
in  Chapter  V,  that  since  this  mesh  does  not  closely  approximate  the  input  objects  per¬ 
imeter  the  resulting  Green's  function  contour  integrals  and  subsequent  far  field  calcu¬ 
lations  have  reduced  accuracy. 

The  nodes  of  the  nodal  rows  form  the  vertices  of  triangular  elements.  An  ele¬ 
ment,  as  defined  in  Chapter  II,  has  three  vertices,  each  assigned  a  unique  number  (1,  2 
or  3).  The  ordering  of  these  vertices  depends  on  the  mesh  orientation  attribute.  Each 
element  also  has  an  associated  relative  dielectric  constant  (e,)  and  relative  permeability 
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ROW  1 


ROW  3 

Figure  14.  Three  Row  Cylinder 

6.  Sorter 

This  subroutine  generates  a  complete  mesh  row  and  all  the  element/node  inter¬ 
connection  relationships.  With  the  nodal  structure  in  place,  individual  nodes  can  be 
assigned  to  individual  elements  within  the  global  mesh  structure.  Each  element  in  a 
given  row  is  assigned  a  unique  local  element  number  starting  with  the  left  most  element 
(relative  to  the  bisection  segment).  See  Figure  15. 


Figure  15.  Element  Numbering  Scheme 


It  is  vital  that  a  method  of  determining  which  nodes  form  the  vertices  of  a  given  element 
and  which  elements  have  a  vertex  attached  to  a  given  local  node.  The  nodes  that  are 
not  part  of  the  offset  contour  have  unknown  field  values.  A  single  node  can  be  con¬ 
nected  to  as  many  as  six  elements,  only  four  of  which  can  be  in  the  current  row.  1  his 
relationship  is  shown  in  Figure  16(a).  This  hexagonal  arrangement  of  six  elements  is 
very  common  within  the  mesh.  Along  the  bisection  segment  or  where  the  mesh  orien¬ 
tation  attribute  changes  from  one  row  to  another,  this  pattern  is  disrupted.  An  example 
of  an  extreme  case  (the  center  of  a  circular  cylinder)  is  shown  in  Figure  16(b).  After  all 
elements  and  nodes  in  an  elemental  row  arc  assigned,  an  ordered  sweep  is  conducted  of 
that  elemental  row  to  determine  which  nodes  arc  connected  to  which  elements,  which 
elements  are  connected  to  which  nodes,  and  how  many  elements  a  single  node  is  con¬ 
nected  too.  A  complete  row  has  now  been  generated.  See  Figure  1 7  on  page  34. 


ATTACHED  TO  ATTACHED  TO 

THIS  NODE  THIS  NODE 

a.  b. 

Figure  16.  Two  Possible  Element  Intersections,  (a)  extreme,  (b)  normal 

The  information  necessary  to  generate  the  rows  and  elements  will  be  used  again,  in 
Chapter  IV,  to  solve  the  (FEBV)  problem. 

7.  Finder 

This  subroutine  determines  the  x,  y  coordinates  of  each  node.  This  data  is  also 
necessary  to  solve  the  FEBV  program  of  Chapter  IV  and  provides  the  plotting 
coordinates  for  the  PLT.DAT  file.  This  process  is  repeated  for  all  rows.  Thus,  the  entire 
object  is  generated  as  the  compilation  of  elemental  rows  made  of  triangles.  See 
Figure  18  on  page  34. 
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Figure  18.  Global  Mesh  Structure 

A  file  is  now  available  for  display  using  a  commercially  available  program  called 
"CURVE-DIGITIZER".  Any  program  that  can  accept  x,  y  coordinate  data  will 
accomplish  the  same  display  process.  Minor  changes  to  the  MESH  program  provided 
in  Appendix  B  may  be  necessary  since  several  "CURVE-DIGITIZER"  plotting  codes 
are  embedded  in  the  file  generation  code. 
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D.  OPTIMIZATION  OF  THE  MESH 

For  most  objects,  it  is  recommended  that  the  MESH  program  be  executed  in  the 
mesh  generation  only  mode  (Field  6  set  to  "M"  or  "m")  prior  to  the  execution  of  the  total 
finite  element  program.  This  will  ensure  that  the  desired  mesh  structure  is  obtained  prior 
to  solving  for  the  unknown  field  values.  The  MESH  (generation  only)  program  requires 
only  a  few  seconds  for  even  the  most  dense  mesh  structures.  It  is  readily  seen  that  as 
the  mesh  resolution  increases,  the  number  of  rows  increases  linearly  while  the  number 
of  unknowns  increases  geometrically.  Therefore,  the  mesh  calculation  times  may  not 
change  appreciably  when  the  mesh  is  made  more  dense,  however,  the  calculation  time 
for  the  finite  element  program  will  rise  geometrically.  For  this  and  other  reasons,  the 
mesh  density  should  be  kept  low.  This  will  lead  to  a  smaller  number  of  unknowns  and 
result  in  faster  program  execution  times.  Figure  19  on  page  36  plots  this  relationship 
for  a  circular  cylinder.  The  solution  for  these  unknown  field  values  will  be  a  topic  of 
Chapter  IV. 

Six  difTerent  mesh  generation  methods  are  available  to  create  mesh  structures.  Some 
of  these  methods  were  evolutionary  in  nature  and  provide  limited  practical  benefit. 
Methods  1  and  6  are  by  far  the  most  useful.  Each  method  is  discussed  below. 

Method  1  constructs  the  bisection  segment  by  connecting  the  first  input  data  point 
to  the  midpoint  of  the  perimeter.  This  segment  is  divided  into  equal  length  segments 
each  separated  by  a  node.  This  method  is  useful  for  very’  simple  objects  such  as  the 
circular  cylinder  shown  in  Figure  20  on  page  37.  This  circular  cylinder  will  be  used  in 
the  explanations  of  all  mesh  generation  methods.  These  illustrations  are  not  intended 
to  optinuze  the  mesh  structure  but  rather  allow  for  easy  comparison  of  the  6  basic 
methods. 

Method  2  constructs  the  bisection  segment  as  described  in  method  1.  The  bisection 
segment  nodes  are,  however,  ordered  difierently.  This  segment  is  divided  by  connecting 
line  segments  from  corresponding  nodes  on  the  left  and  right  side  perimeter.  Where 
these  line  segments  intersect  the  bisection  segment,  a  node  is  placed.  This  leads  to  un¬ 
equally  spaced  bisection  segments.  This  can  be  seen  in  Figure  21  on  page  37. 
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Figure  19.  Uiiknouns  and  Number  of  Rons  Versus  Mesh  Resolution 


Method  3  modilies  method  1  by  allowing  the  user  to  specify,  using  Field  17,  the  stop 
node  for  the  bisection  segment.  This  method  can  be  useful  to  rapidly  adjust  around 
slightly  irregular  objects.  This  can  be  seen  in  Figure  22  on  page  38. 

Method  4  combines  methods  2  and  3  by  using  the  connected  line  segment  technique 
of  method  2  to  determine  the  node  positions  on  the  bisection  segment,  the  stop  node  of 
which  is  specified  as  in  method  3.  This  can  be  seen  in  Figure  23  on  page  39. 

Method  5  improves  upon  method  4  by  repositioning  the  unequally  spaced  bisection 
nodes.  Linearly  interpolated  positions  for  the  nodes  leads  to  equal  spacing.  This 
method  reduces  the  node  "bunching"  that  frequently  occurs  with  methods  3  and  4.  This 
can  be  seen  in  Figure  24  on  page  40. 
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Figure  20.  Method  1  Mesh  Structure  Example 
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Figure  21.  Method  2  Mesh  Structure  Example 


Figure  22.  Method  3  Mesh  Structure  Example 


Method  6  is  the  final  improvement  in  which  method  5  was  modified  to  allow  for  a 
user-specified  start  node.  This  provides  the  user  with  the  ability  to  start  and  stop  the 
bisection  segment  at  any  input  node  without  having  to  rearrange  the  input  data.  Since 
method  6  contains  all  of  the  capabilities  of  the  other  five  methods,  it  is  almost  exclu¬ 
sively  used  for  mesh  generation.  This  can  be  seen  in  Figure  25  on  page  40.  There  are 
situations  that  could  be  best  served  by  one  of  the  other  methods.  For  example,  a  cir¬ 
cular  cylinder  and  other  simple  symmetric  objects  can  be  represented  using  method  1. 
Method  2  would  work  well  w'ith  a  square  or  diamond  shape.  These  objects  w'ould  be 
bisected  by  a  diagonal.  Method  3  would  be  most  suited  for  a  symmetric  object  with  a 
planar  material  interface.  The  more  dense  mesh  would  be  used  for  the  higher 
permittivity  material.  Method  4  would  work  well  with  a  rhombus  or  other  slightly 
asymmetric  object.  Method  5  is  almost  as  versatile  as  method  6,  and  W’ould  work  well 
in  any  situation  where  the  start  node  is  fixed.  A  great  deal  of  planning  is  not  needed  in 
designing  most  mesh  structures  since  they  are  calculated  and  displayed  in  a  matter  of 
seconds. 
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Iterative  selection  of  different  mesh  generation  parameters  has  proven  to  be  the  best 
technique.  A  summary  of  all  mesh  generation  capabilities  is  provided  in  Table  1. 

Appendix  C  contains  a  program  called  READ.FOR.  This  program  takes  the  output 
data  file  from  the  "Curve-Digitizer"  CAD  program,  called  F1NALDVVG.DAT,  and  after 
receiving  the  answers  to  several  prompted  questions,  creates  a  new  INPUT.DAT  file. 
Typically,  the  CAD  program  is  used  to  generate  the  perimeter  of  the  object.  This  may 
be  as  a  series  of  points,  line  segments  or  a  combination  of  the  tw’o.  The  answers  to  the 
prompted  questions  provide  the  additional  information  needed  to  fill  the  remaining  data 
fields.  This  is  intended  to  be  a  first  step  towards  allowing  a  user  to  specify  or  design  an 
object  and  then  be  able  to  calculate  the  scattered  fields  from  this  object.  Although  it  is 
far  from  efficient,  it  does  serve  a  definite  purpose.  Further  refinement  of  this  program 
will  allow  for  an  easier  user  interface.  This  should  enable  technically  trained  personnel, 
without  a  detailed  understanding  of  these  programs,  to  benefit  from  the  Field  Feedback 
Formulation. 


Figure  24.  Method  5  Mesh  Structure  Example 


Figure  25.  Method  6  Mesh  Structure  Example 


Table  1.  SUMMARY  OF  MESH  GENERATION  CAPABILITIES 


Method 

Number 

Straight 

Line 

Bisection 

Straight 
Line  From 
Left  to 
Right  Side 

User  Se¬ 
lected  Stop 
Node 

Equally 

Spaced 

Bisection 

Nodes 

User  Se¬ 
lected  Start 
Node 

1 

X 

X 

2 

X 

3 

X 

X 

4 

X 

X 

5 

X 

X 

6 

X 

X 

X 
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IV.  FINITE  ELEMENT  BOUNDARY  VALUE  PROGRAM 


A.  INTRODUCTION 

The  Finite  Element  Boundan.'  Value  (FEBV)  program  is  the  feed  forw'ard  (U)  oper¬ 
ator  in  the  Field  Feedback  Formulation  (F^)  as  shown  in  Figure  2  on  page  4.  This 
program  solves  the  Helmholtz  equation,  as  discussed  in  Chapter  II,  for  the  Dirichlet 
boundary  condition  specified  in  the  input  data  file,  as  discussed  in  Chapter  III.  These 
boundary  conditions  are  imposed  on  the  offset  boundary^  contour.  The  purpose  of  this 
program  is  to  find  the  unknown  field  values  inside  the  offset  boundary  contour  within 
the  input  object.  Since  the  ultimate  goal  is  to  obtain  the  scattered  far  fields  from  the 
object,  the  unknown  field  values  on  the  perimeter  are  of  primary  interest.  It  will  also 
be  necessary  to  approximate  the  normal  derivative  of  the  field  at  the  object  perimeter. 
As  derived  in  Chapter  II,  the  goal  of  this  program  is  to  solve, 

-H  =  -  [?,]'F, 

The  A  matrix  represents  the  effect  that  the  I  —  1  row  has  on  the  /'*  row  values.  Similarly, 
the  B  matrix  represents  the  effect  that  the  I-th  row  has  on  itself,  while  the  C  matrix  re¬ 
presents  the  efiect  that  the  I  +  1  row  has  on  the  I-th  row.  All  other  rows  do  not  effect 
the  I-th  row.  This  is  due  to  the  pyramidal  basis  functions,  discussed  in  Chapter  II,  that 
all  have  zero  value  when  a  node  is  not  directly  connected  to  another  node.  The  P  matrix 
represents  the  combined  effects  of  the  boundary  nodes  on  the  1-th  row.  These  values 
are  transposed  to  the  right  side  of  the  equality  to  form  the  system  forcing  function.  It 
is  worth  remembering  that  for  the  I  =  I  row,  the  A  matrix  =  0  and  that  for  the  I  = 
I. MX  row,  the  C  matrix  =  0.  Thus,  using  the  row  by  row  stepping  process,  first  dis¬ 
cussed  in  Chapter  III  of  the  mesh  generation  process,  the  A,  B,  C  and  P  matrices  can 
be  filled.  There  is  an  A,  B,  C  and  P  matrix  for  each  row,  and  it  is  necessary'  to  calculate 
the  functionals  for  two  elemental  rows  to  fill  one  set  of  matrices.  It  is,  therefore,  obvious 
that  the  data  is  used  twice,  once  for  the  current  row  and  again  for  the  next  row.  Spe¬ 
cifically,  the  functionals  necessary  to  complete  the  B  matrix  and  totally  fill  the  C  matrix 
is  used  again  to  totally  fill  the  next  row's  A  matrix  and  partially  fill  the  B  matrix.  Thus, 
for  a  row  I  -  1,  (such  that  I  -  1  1),  all  of  the  A  matrix  and  part  of  the  B  and  P  matrices 

are  filled.  When  the  row  is  stepped  to  I,  the  B  and  P  matrices  are  completed  and  all  of 
the  C  matrix  is  filled.  The  A,  B,  C  and  P  matrices  represent  tri-block  matrices  within  a 


much  larger  global  matrix.  This  row  of  tri-block  matrices  in  the  global  matrix  can  be 
partially  solved  using  the  forward  portion  of  the  Ricatti  transform.  The  Ricatti  trans¬ 
form  is  a  numerical  technique  that  optimizes  the  solution  for  banded  (tri-block)  systems 
of  linear  equations.  The  equations  are, 

^;+i  =  —  {.Bi  A-  AiR^  *  Q 


and 


=  {B,  =  (5;  AiRf'  •  (P  -  AiS,) 

where  is  the  i'*  -I-  I  R  matrix  and  the  S,.,  is  the  /'*  -I-  1  S  vector.  Note  that  the  next 
rows  R  matrix  and  S  vector  depends  on  the  previous  rows  R  matrix  and  S  vector. 
During  the  forward  step  (MARCH  subroutine)  an  R  matrix  and  S  vector  is  generated 
for  each  row.  When  all  rows  have  been  calculated,  a  back  sweep  computes  the  unknown 
field  values,  The  equation  is 


^t-x  =  Ri^i  +  Sf. 

This  back  sweep  must  read  the  RS  data  from  the  disk  backwards.  This  is  a  storage  in¬ 
tensive  process  for  all  but  the  most  trivial  problem.  The  field  values  is  first  found 
by  remembering  that  the  last  row  (  I  =  I. MX  ),  C,.„^  =  0.  With  >  Bnix-\  = 

and  Sjxix-]  .  which  is  the  last  calculated  S  vector.  The  recursion  continues  until 

all  of  the  ij/i 's  have  been  calculated. 

B.  FINITE  ELEMENT  BOUNDARY  VALUE  PROGRAM 

The  eight  subroutines  comprising  the  FEBV  program  are  discussed  below. 

1.  Zero 

This  subroutine  fills  all  the  array  positions  of  the  A,  B,  C  and  P  matrices  wth 
0  +  jO.  This  zeroing  is  necessar>'  after  all  calculations  concerning  a  given  row  are 
completed. 

2.  Varint  (Variational  Integration) 

This  subroutine  calculates  the  complex  functionals  for  a  given  input  element. 
The  functionals  reflect  the  effect  that  each  node  has  on  the  three  nodes  associated  with 
an  element.  The  subroutine  must  be  provided  with  the  X,  Y  coordinates  of  the  element 
nodes,  and  the  material  parameters  e,  and  n,.  These  functionals  are  returned  in  a  3  x  3 
complex  matrix.  The  area  of  the  input  element  is  also  provided  as  a  by-product  of  this 
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numerical  integration.  The  areas  of  each  element  are  sorted  and  used  to  find  the  largest 
and  smallest  elements.  Once  found,  the  largest  and  smallest  areas  are  combined  to  form 
the  area  ratio,  which  is  defined  as. 


area  ratio  = 


maximum  area 
minimum  area 


This  ratio  should  be  kept  as  small  as  possible.  For  a  circular  cylinder,  values  slightly  less 
than  2.0  are  possible.  For  more  complicated  objects,  the  ratio  can  get  considerably 
larger.  A  ratio  ^  2.5  causes  a  screen  message  of  "YOU  SHOULD  CONSIDER 
ABORTING  THIS  RUN  AND  LOOKING  AT  THE  MESH  IN  CURVE  DIGITIZER. 
A  BETTER  METHOD  MAY  BE  .AVAILABLE".  It  may  or  may  not  be  possible  to 
obtain  an  area  ratio  <  2.5.  The  intention  of  the  area  ratio  is  to  insure  that  a  uniform 
mesh  is  constructed  prior  to  attempting  a  problem  solution.  Smaller  area  ratios  are  in¬ 
dicative  of  mesh  structures  that  do  not  have  grossly  different  element  sizes.  This  will 
lead  to  more  accurate  finite  element  solution. 

3.  Fill 

This  subroutine  calculates  and  stores  all  of  the  functionals  for  an  entire  ele¬ 
mental  row.  This  is  accomplished  by  repeated  VARINT  calls  for  each  element  that 
comprises  an  elemental  row. 

4.  BNDC  (Boundary’  Condition) 

This  subroutine  calculates  and  stores  the  boundary^  conditions  desired  for  the 
offset  boundary’  contour.  These  conditions  can  be  plane  wave,  cylindrical  modes  or  in¬ 
put  manually.  For  plane  wave  conditions,  the  percent  error  of  the  FEBV  program  will 
also  be  returned.  This  calculation  only  has  meaning  for  the  uniform  case  (Field  5  set  to 
"U"  or  "u").  .Memory’  limitations  do  not  allow’  this  feature  for  the  cylindrical  mode 
boundary  conditions  for  modes  other  than  zero  and  one.  A  separate  routine  could  be 
made  available  for  offline  percent  error  calculations.  No  such  feature  is  provided  if  the 
boundary  conditions  are  manually  input.  Any  desired  boundary  condition  may  be  gen¬ 
erated  by  modifying  or  appending  the  proper  code  to  this  subroutine. 

5.  Loader 

This  subroutine  loads  the  A,  B,  C  and  P  matrices  for  a  given  row.  For  all  but 
the  I  =  1  and  I.MX  rows,  tw’O  calls  of  FILL, 'VARINT  for  tw’o  elemental  row's  of  data 
are  required  to  fill  the  A,  B,  C  and  P  matrices.  This  is  an  additive  process  that  starts 
after  the  ZERO  subroutine  initializes  the  matrices.  Successive  functionals  are  added  to 
the  existing  data  in  the  proper  array  positions. 
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6.  March 

This  subroutine  performs  the  forward  portion  of  the  Riccati  transform.  The 
output  data,  called  RS.DAT  (R  matrix  and  S  vector),  are  stored  on  disk. 

7.  CSMINV  (Complex  Square  Matrix  Inversion) 

This  subroutine  accomplishes  the  complex  matrix  inversion  required  by  the 
forw’ard  Riccati  transform.  A  maximum  dimension  of  50  x  50  was  established  to  limit 
memory  utilization.  This  allows  for  a  maximum  unknown  width  of  50. 

8.  Sweep 

All  of  the  above  subroutines  are  called  at  least  once  for  each  row'.  The  sw'eep 
routine  is  called  only  after  all  of  the  row  calculations  are  completed,  and  then  only  once. 
This  subroutine  conducts  the  Riccati  back  sw'eep  by  reading  the  RS  data  generated  by 
the  MARCH  subroutine.  This  data  is  read  off  the  disk  backwards,  using  the 
FORTRAN  "backspace"  command.  The  returned  field  values  are  actually  individual 
contributions  due  to  a  unit  valued  basis  function  being  individually  applied  to  each  of 
the  boundar}'  nodes.  In  so  doing,  the  problem  need  only  be  solved  once.  After  the  data 
is  stored,  in  matrix  form  as  the  U.DAT  file,  the  boundary  value  problem  may  be  solved 
for  any  incident  field  by  multiplying  the  U  matrix  by  the  new'  incident  fields.  Thus, 

perimeter  incident- 

A  summary  of  the  input,  output  and  error  data  is  provided  in  the  OUTPUT.DAT  file. 

9.  Save 

This  subroutine  was  necessary  to  allow'  for  reasonably  sized  problems.  Since 
all  of  the  Field  Feedback  Formulation  code  could  not  fit  into  640  kilobytes  of  memory, 
the  program  was  divided  into  two  parts.  All  of  the  data  necessary'  to  perform  the  pro¬ 
grams  is  saved  in  the  F3.DAT  file.  This  data  is  then  read  by  the  field  feedback  program. 
This  technique,  though  very  inefficient,  circumvents  the  640  kilobyte  memory  limitation 
imposed  by  the  IBM  Disk  Operating  System  (DOS).  Efibrts  to  convert  this  code  to  run 
under  a  compiler  that  does  not  have  a  640  kilobyte  memory  limitation,  such  as  Micro¬ 
way  NDP  FORTRAN  compiler,  w'ill  be  pursued  at  a  later  date. 

C.  VARINT  VALIDATION 

The  first  step  in  the  program  validation  required  an  understanding  of  the  error 
convergence  of  the  variational  elements  as  a  function  of  element  size,  d  and  material 
properties,  c,  and  ,u.  •  The  test  program  provided  in  Appendix  D  varies  the  element  size, 
in  wavelengths,  for  a  test  mesh  structure.  This  test  structure  shown  in  Figure  26  on 
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page  46  has  only  one  unknown  at  the  center.  The  mesh  is  made  of  four  adjacent  ele¬ 
ments.  These  elements  are  inside  a  uniform  material  having  properties,  e,  and  1  he 
other  four  nodes  are  established  as  boundary  nodes. 


PLANEWAVE 


Figure  26.  Test  Mesh  Structure 

A  plane  wave,  of  user  specified  frequency  and  amplitude,  is  used  to  determine  the 
boundarj'  conditions  on  the  four  boundary  nodes.  This  plane  wave  is  propagating  down 
the  vertically  oriented  axis.  The  unknown  nodal  value  is  calculated  and  compared  to  the 
actual  plane  wave  value  ( I  +  jO).  A  percent  error  is  calculated  as  the  dimension  of  the 
elements,  d  are  reduced.  As  expected,  the  solution  became  more  accurate  as  the  ele¬ 
ments  become  smaller.  A  plot  of  the  results  of  this  test  program  is  provided  in 
figure  27  on  page  47.  This  data  was  generated  with  e,=  I  -i- JO,  and  fx,  —  I  +  JO. 
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Figure  27.  Solution  Error  for  a  Test  Mesh  Structure 

The  region  of  interest  is  where  the  error  approaches  zero.  An  expanded  plot  of  this  re¬ 
gion  is  provided  in  Figure  28  on  page  48.  The  accuracy  of  the  solutions  were  desired 
within  1  percent  of  the  exact  value.  I  his  requires  an  element  that  is  ^  ,  in  tlie  ma¬ 

terial.  To  ensure  the  desired  error  was  achieved,  elements  were  usually  scaled  to  be 
S  ,  in  the  material.  The  phase  of  the  plane  wave  was  also  varied  to  determine  the 
elfect  on  convergence.  No  significant  effect  was  noted. 


Figure  28.  Solution  Error  for  a  Test  Mesh  Structure  (Expanded) 


D.  FINITE  ELEMENT  BOUNDARY  VALUE  PROGRAM  VALIDATION 

1  he  next  validation  step  required  an  actual  object  to  calculate  the  fields  in  and  on. 
The  simplest  object  ;vas  actually  not  an  object  at  all,  but  rather  an  imaginary  circular 
cylinder  in  free  space.  A  plane  wave  was  propagated  through  this  free  space.  Since  no 
material  interface  existed,  the  exact  solution  would  be  known  for  all  positions.  The 
plane  wave  established  the  boundary  conditions  on  the  offset  contour.  Trial  runs  were 
conducted  varying  the  mesh  resolution  and  boundary  offset  contour  distance.  With  the 
error  defined  as, 


error  = 


Ifcalculated  value  —  actual  value)^ 


Z(actual  value)^ 


two  errors  were  defined.  The  perimeter  error  only  considers  the  perimeter  nodes  in  the 
error  calculation.  The  bisection  segment  error  only  considers  the  bisection  segment 
nodes,  with  the  exception  of  the  two  end  nodes,  in  the  error  calculation.  The  end  nodes 
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of  the  bisection  segment  are  part  of  the  perimeter  and  are,  therefore,  not  considered. 
As  expected,  the  perimeter  error  could  be  rapidly  reduced  by  decreasing  the  offset  dis¬ 
tance.  This  is  due  to  the  fact  that  the  node  or  nodes  closest  to  an  unknown  node  dom¬ 
inate  the  field  contributions  at  this  node.  Again,  a  goal  of  <  1%  error  was  desired.  The 
reduced  offset  distance  had  a  very  small  effect  on  the  bisection  segment  error.  The  only 
way  to  significantly  reduce  this  error  was  to  increase  the  mesh  resolution.  An  increased 
mesh  resolution  reduced  both  errors.  Figure  29  show's  the  perimeter  error  as  a  function 
of  the  number  of  bisection  segment  nodes. 


Figure  29.  Perimeter  Error 

The  three  curves  are  for  offset  distances  of  0.05,  0.025  and  0.01  X  .  Figure  30  on  page 
50  show's  the  bisection  segment  error  as  a  function  of  tlie  number  of  bisection  segment 
nodes.  The  two  curves  are  for  offset  distances  of  0.05  and  0.025  X^ .  For  offset  distances 
less  than  0.025  2o ,  the  curves  are  almost  identical  to  the  0.025  curve.  These  curves 
are  omitted  for  clarity. 
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Figure  30.  Bisection  Segment  Error 


The  calculated  and  exact  field  values  for  a  0.5  diameter  circular  cylinder,  with  a 
O.O5.I0  mesh  resolution  and  oflset  distance  and  c,  =  1  +/)  is  shown  in  Figure  31  on  page 
51.  The  exact  fields  values  for  the  real  and  imaginary  portion  of  the  plane  wave  are 
shown  as  squares  and  diamonds,  respectively.  The  solid  curves  are  the  calculated  field 
values.  The  perimf  r^r  and  bisection  errors  were  0.74  and  1.69  percent,  respectively. 
Figure  32  on  page  52  shows  the  efiects  of  not  properly  adjusting  the  mesh  resolution 
and  olfset  distance  when  the  material  is  changed.  In  this  case,  the  permittivity  was 
changed  to  c,  =  4  +  jO.  This  caused  the  perimeter  and  bisection  segment  errors  to  in¬ 
crease  to  18.2  and  41.1  percent,  respectively.  The  mesh  resolution  and  offset  distance 
were  reduced  to  0.0212o  and  the  permittivity  was  changed  to  c,  =  4  —  jA.  The  results  are 
shown  in  Figure  33  on  page  53  .  This  proper  selection  of  mesh  resolution  and  offset 
distance  reduced  the  perimeter  error  and  bisection  segment  error  to  0.44  and  0.56  per¬ 
cent,  respectively.  Note  that  in  the  lossy  material  case  the  two  errors  are  very  close  in 
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magnitude.  This  is  due  to  the  way  that  the  error  is  defined.  This  definition,  in  a  lossy- 
material,  overemphasizes  the  objects  leading  edge. 


fwlitwNr  NumNr 

Figure  31.  0,5  ).  cylinder,  e,  =  1  +  jO 

This  portion  of  the  object  (as  well  as  the  trailing  edge)  is  the  most  accurate  for  the 
bisection  segment  calculations.  This  is  because  of  the  relative  closeness  of  these  nodes 
to  the  perimeter.  For  the  lossless  material,  the  bisection  segment  error  is  typically  two 
to  three  times  the  perimeter  error,  for  the  same  number  of  bisection  nodes.  The  con¬ 
clusion  that  the  offset  distance  should  be  made  as  small  as  possible  in  order  to  minimize 
the  perimeter  error  is  correct,  but  will  prove  to  be  counterproductive  in  the  long  run. 
There  is  a  competing  effect  that  will  require  this  contour  offset  distance  to  be  as  large 
as  possible.  This  effect,  associated  with  the  accuracy  of  the  Green's  function  contour 
integral,  will  be  discussed  in  Chapter  V. 
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Figure  33.  0.5  X  cylinder,  t,  -  4  —  \4 


E.  INHOMOGENEITY 

Until  this  point,  all  testing  was  done  in  circular  homogeneous  dielectric  materials. 
It  was,  at  a  minimum,  desirable  to  place  the  object  in  a  vacuum  to  calculate  the  scattered 
fields.  This  requires  the  first  two  and  last  two  elements  of  each  row  to  have 
c,  =  I  +J0  and  ti,  =  1  +7O.  These  elements  represent  the  space  between  the  objects  per¬ 
imeter  and  the  offset  boundary  contour.  Since  the  plane  wave  solution  is  no  longer  valid 
for  this  geometry,  a  new  problem  with  a  known  solution  was  needed. 

Given  a  homogeneous  dielectric  circular  cylinder  of  radius  a  and  permittivity  t, ,  as 
shown  in  Figure  34  on  page  54,  it  can  be  shown  that  [Ref.  8  ], 

4>)  =  A„J„{k,R)  cos  n<i>. 
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where, 


and 


^„{R,  0)  is  the  field  value  inside  the  dielectric  material. 

^  ^  Ji^R,)[j„(k,RM^{R,)  -  krJ’„{WH^^\Rj]  - 
"  Hj^\R,)lJ„{WJ’„iR„)  ~  kr/'„iWJM2 

J„  ==  Bessel  Function  of  the  l”  kind  of  order  n 

K  ~ 
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An  additional  formulation  is  required  for  the  fields  outside  of  the  cylinder.  A  cylindrical 
wave  (mode,  n  =  1)  was  applied  to  a  homogeneous  dielectric  cylinder.  The  exact  and 
finite  element  field  values  were  calculated  and  compared.  Figure  35  on  page  55  ,  shows 
a  typical  result.  Mesh  resolution,  offset  distance  and  material  permittivity  were  varied 
to  determine  convergence.  The  results  were  similar  to  the  dielectric  homogeneous  plane 
wave  case  discussed  earlier. 


Parimalw  Nods  NiimlMr 

Figure  35.  Typical  Cylindrical  Mode  Boundary  Value  Problem  Result 


F.  FINITE  ELEMENT  CONCLUSIONS 


High  accuracy  solutions  can  be  obtained  by  maintaining  a  mesh  resolution  of 
<~^  in  the  material.  The  offset  distance  should  be  kept  at  approximately  the  same 
magnitude  as  the  mesh  resolution,  but  may  be  as  large  as  .  Increased  accuracy  re¬ 


quires  an  increased  mesh  resolution.  This  increase  in  accuracy  is  at  the  expense  of  in¬ 
creased  processor  time. 
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V.  FIELD  FEEDBACK  PROGRAM 


A.  INTRODUCTION 

The  Field  Feedback  Program  is  the  feedback  (T)  operator  in  the  f\  This  program 
utilizes  the  output  of  the  finite  element  program  and  the  input  incident  fields  to  evaluate 
a  "near  field"  Green's  function  integral.  This  integral  is  called  "near  field"  in  that  the 
integral  will  be  performed  within  of  the  object  perimeter.  As  seen  in  Figure  36, 
three  surface  contours  are  defined.  The  boundary  contour  is  where  the  incident  field 
boundary  condition  is  applied. 


Figure  36.  Contour  Arrangement 


The  perimeter  contour  is  where  the  finite  element  program  solves  for  the  field  values  on 
the  object  perimeter.  The  geometric  contour  is  midway  between  the  boundary  and  the 
perimeter  and  is  the  contour  over  which  the  Green's  function  contour  integral  is  per- 
formed.  This  is  necessary  to  allow  the  to  be  approximated  by  a  finite  difl'erence 
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technique.  The  ^cc  is  the  average  of  the  field  values  on  the  boundary'  and  perimeter. 
Thus, 


d\l/  boundary  ^perimeter) 

6n  8n 


and 


^cc  — 


perimeter  "b  ^boundary) 


The  Field  Feedback  program  is  provided  as  Appendix  E. 

B.  FIELD  FEEDBACK  PROGRAM 

1.  Input 

This  subroutine  reads  the  finite  element  data  stored  in  the  F3.DAT  file  by  the 
SAVE  subroutine.  All  necessaiy^  nodal  X,  Y  coordinates  are  calculated. 

2.  TMAT  (T  Matrix) 

This  subroutine  loads  a  complex  matrix  called  TMAT.  Each  element  of  the  T 
matrix  is  the  result  of  a  Green's  function  contour  integral.  This  integral  is  conducted 
on  the  GC  contour.  The  m-th  column  of  the  T  matrix  is  evaluated  with  a  single  unit 
valued  basis  function  on  the  m-th  boundarj-  node.  The  equation, 


^n,m  = 


f 

r 

G-rr^ 

on 

,  dG  ~ 

-  ^n^~ 

on 

'gc 

dOC, 


may  now  be  numerically  evaluated  at  each  of  the  n  boundary  nodes.  This  process  is 
repeated  for  each  row  until  the  entire  matrix  is  filled. 

The  numeric  integration  uses  a  rectangular  mid-point  approximation  technique. 
For  this  linearized  problem,  this  is  equivalent  to  a  trapezoidal  integral  technique.  The 
integral  path  between  any  given  two  nodes  is  subdivided  based  on  the  distance  to  the 
first  point.  This  subdivision  maybe  controlled  by  three  input  variable  fields.  A  multi¬ 
plicative  scale  factor  and  offset  terms  are  available.  To  date,  the  utilization  of  the  scale 
factor  (set  >  1.0)  and  the  offset  terms  (set  >  0  )  have  not  proved  necessaiy.  This  may¬ 
be  necessary  for  more  <"omplicated  geometric  structures. 
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3.  CNSOLV 

This  subroutine  solves  the  equation. 


Qi  ~  '^matrix^  *  ^boundary* 


where  C„  is  the  total  scattered  fields  on  the  boundary,  I  is  the  identity  matrix  and 
^boundcry  is  the  initially  specified  incident  fields. 

4.  FFLD  (Far  Fields) 

This  subroutine  calculates  the  scattered  far  field  or  bistatic  radar  cross  section 
per  unit  length  of  the  object.  Any  desired  integer  angular  resolution  for  the  calculation 
may  be  specified.  The  final  results  are  stored  in  the  FFPAT.DAT  file. 

C.  FIELD  FEEDBACK  VALIDATION 

Given  a  plane  wave  boundary  conditions  imposed  on  an  imaginary  cylinder  in  free 
space,  a  Green's  function  contour  integral  may  be  performed  around  the  cylinder.  This 
cjdinder  is  imaginary  in  the  sense  that  it  does  not  actually  exist.  The  cylinder  is  actually 
an  artificial  boundary  that  does  not  form  a  material  interface.  For  points  inside  this 
cylinder,  the  resulting  integral  should  equal  the  negative  of  the  actual  plane  wave  value 
at  that  point  in  freespace.  For  points  outside  the  cylinder,  the  resulting  integral  should 
equal  zero.  Several  test  conditions  were  investigated  to  ensure  that  the  numerical  inte¬ 
gration  of  the  Green's  function  did  arrive  at  the  expected  values.  Maximum  absolute 

errors  for  these  tests  were  less  than  7  x  10"’ .  These  test  cases  started  with  exact  values 

dll/ 

of  the  field,  if/  and  the  normal  derivative,  .  After  initial  validation,  finite  element 
dij/ 

calculated  il/  and  values  were  used.  Errors  increased  by  approximately  a  factor  of 
10.  Numerous  competing  effects  were  observed  with  two  dominant  effects  becoming 
evident. 

1.  Small  Object  Phenomenon 

For  very  small  objects,  where  only  a  few  elements  are  needed  to  meet  the  -~ 
mesh  resolution  requirement,  the  normal  derivative  accuracy  is  crucial.  Since  the  normal 
derivative  is  the  calculated  normal  to  the  piecewise  linear  approximation  of  the  objects 
perimeter,  this  fit  is  usually  not  adequate.  To  prevent  this  problem,  additional 
perimeter/boundary  nodes  are  needed.  This  is  easily  accomplished  by  increasing  the 
mesh  density,  which  is  controlled  by  the  mesh  resolution  (input  field  7).  For  these  very 
small  objects,  the  RCS  is  almost  uniform  for  the  T.M  case  and  co-sinusoidal  for  the  TE 
case.  The  utility  of  this  calculation  must,  therefore,  be  questioned. 
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2.  Offset  Distance  Phenomenon 

The  offset  distance  must  not  be  made  too  small  or  the  accuracy  of  the  Green's 
function  integral  will  suffer.  For  a  circular  cylinder,  the  optimal  offset  was  between  0.03 
A  and  0.035  /.  Note  that  these  distances  are  always  <  . 
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VI.  VALIDATION 


A.  INTRODUCTION 

The  difficulty  in  the  total  Field  Feedback  Formulation  (f ’)  program  validation  is  in 
finding  problems  that  have  well  established  or  accepted  solutions.  The  total  program 
is  the  combination  of  the  mesh  generation/finite  element  program  and  the  field  feedback 
program.  If  possible,  a  problem  that  has  a  closed  form  analytical  solution  is  desired. 


B.  HOMOGENEOUS  CIRCULAR  CYLINDRICAL  SCATTERING 

This  is  the  first  test  case  for  the  program.  A  penetrable  circular  cylinder  has  a 
closed  form  analytic  solution,  so  the  final  results  could  be  verified  against  exact  values. 
It  can  be  shown  that  for  the  TM  case  (E  -  wave)  the  reflection  coefficient  is, 


and  that  for  the  TE  case  (H  -  wave)  the  refection  coefficient  is, 

UWJ'M  -  7-r 

j-TM  _  __L _ ^  ’’ _ J_ 

where  k,  =  ,  y,  =  \  \  ^  1  scattered  field  for  the 

TM  case  is, 


4(F,(f>)  =  -£' 


incident 


oo 

To/^'^F)  +  cos  n4> 


/J=l 


Similarly,  the  scattered  field  for  the  TE  case  is, 


H'ziR,  4>)  =  -H' 


incident 


ro/4'’(F)  +  2 


oo 


r.//; 


■<2), 


\R)  cos  rifh 
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In  the  far  field  region,  as  R  approaches  infinity,  the  scattering  width  may  be  defined  as. 


,  ..  Iw’l’ 


and 


o(0)  =  -^  1  To  +  2^r„  cos  n4>  1  ^ 

The  source  code,  without  the  Bessel  function  routines,  for  these  calculations  is  provided 
as  Appendix  F. 

Three  different  radius  dielectric  cylinders  were  tested.  Wavenumber  normalized  radii 
of  0.5,  1.0  and  2.0  with  e,  =  2.56  +  jO  were  selected.  The  TE  and  TM  results  of  the  0.5, 
1.0  and  2.0  radii  problems  are  provided  as  Figure  37  on  page  62,  Figure  38  on  page 
63,  and  Figure  39  on  page  64  .  To  validate  the  Chapter  VI  conclusion  that  the  unit 
normal  was  the  cause  of  the  errors  for  very  small  circular  cylinders,  the  c,  was  increased 
from  2.56  to  25.6.  This  value  still  meets  the  requirement  for  mesh  resolution  not  to  ex¬ 
ceed  .  The  TE  and  T.M  results  are  provided  as  Figure  40  on  page  65  .  Actual  cross 
section  values  are  indicated  by  the  squares  and  the  Field  Feedback  Formulation  sol¬ 
utions  are  plotted  as  a  single  solid  curve. 

C.  HOMOGENEOUS  IRREGULAR  OBJECTS 

The  results  detailed  in  [Ref.  10]  were  next  validated  with  the  program.  This  ob¬ 
ject,  as  seen  in  Figure  41  on  page  66  ,  is  a  dielectric  shell  with  inner  radius  =  .25 
outer  radius  =  .30  Zq  and  s,  =  4  +  jO.  The  comparison  of  the  [Ref  10]  data  and  the 
F  results  are  provided  as  Figure  42  on  page  67. 
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Pill  (In  4«trMt) 

Figure  38.  Cylinder,  TE  and  TM  Case,  k^a  =  1.0,  c,  =  2.56 
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PM  (In  4«arM«) 


Figure  39.  Cylinder,  TE  and  TM  Case,  =  2.0,  8,  =  2.56 


■tiiais 


Ptil  (In 

Figure  40.  Cylinder,  TE  and  TM  Case,  k,fl  -  0.5,  e,  =  25.6 
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PLANEWAVE 


Figure  41.  Dielectric  Shell  Mesh 

The  TE  case  shows  excellent  agreement.  The  TM  case  tends  to  diverge  at  the  smaller 
values  of  0. 

D.  AN  INHOMOGENEOUS  OBJECT 

The  results  detailed  in  (Ref.  11]  were  next  validated  with  the  f’ program.  This  ob¬ 
ject,  as  seen  in  Figure  43  on  page  68,  is  a  dielectric  ring  with  inner  radius  =  .25  X„ ,  outer 
radius  =  .30  J.#  and  c,  =  A  +  jO.  The  exact  solution  to  this  problem  is  available.  [Ref. 
11] 
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Figure  43.  Dielectric  Ring 

The  mesh  generation  program  was  not  modified  to  conform  to  both  the  outer  and  inner 
material  boundaries.  The  original  mesh  generation  program  design  concept  was  to 
closely  fit  the  mesh  to  the  objects  outer  perimeter  and  then  allow  for  material  inhomo¬ 
geneities  to  be  accounted  for  by  diflerent  material  parameters  being  assigned  to  individ¬ 
ual  elements.  A  section  of  the  generated  mesh  is  show  n  in  Figure  44  on  page  69.  Note 
the  additional  curved  line  showing  the  .25  inner  radius.  This  curve  does  not  follow' 
the  established  element  boundaries.  The  effective  ring  is  shown  in  Figure  45  on  page 
69.  This  resulted  in  the  inner  radius  having  an  irregular  pattern  as  elements  near  the 
material  interface  varied  in  material  composition.  This  is  similar  to  the  granular  noise 
problem  characteristic  of  delta  modulation  communication  channels.  The  TE  and  TM 
results  are  provided  as  Figure  46  on  page  70. 
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Figure  44. 


Figure  45. 


Partial  Mesh  with  Inner  Radius  Curve 


Effective  Geometry  for  the  Dielectric  Mesh 
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□  -  TU  Exact 
—  -  I  E  Calculated 


Phi  (bi  dagrtM) 


□  -  TiVl  Calculated 
—  -  TM  Exact 


PM  (In  dagrMi) 


VII.  CONCLUSIONS 


A.  RESULTS 

The  Field  Feedback  Formulation  proved  to  be  an  excellent  tool  for  calculating  the 
internal  and  scattered  fields  from  the  tested  inhomogeneous  asymmetric  objects.  The 
keys  to  satisfactory  results  are, 

•  Keep  the  maximum  element  dimension  , 

•  Maintain  a  mesh  with  a  uniform  distribution,  and 

•  Maintain  the  offset  distance  near  the  mesh  resolution  in  magnitude,  but  no  greater 
than  . 

The  techniques  used  proved  to  be  capable  of  adapting  to  a  wide  variety  of  situations. 
These  adaptations  did  not  require  program  modifications  or  reprogramming.  The  nu¬ 
merous  built-in  features,  as  discussed  in  the  input  field  description  of  Chapter  III,  al¬ 
lowed  for  this  robustness. 

B.  RECOMMENDATIONS  AND  EXTENSIONS 

With  the  basic  program  validated,  future  testing  and  development  should, 

•  Emphasize  large  object  validation  against  accepted  results. 

•  Stress  inhomogeneity  and  irregularity  in  all  testing. 

•  Optimize  the  T  matrix  element  calculation  by  improving  the  Green's  function 
contour  integral. 

•  Evaluate  the  usefulness  of  the  maximum  distance  feature  (Field  12).  Beyond  this 
maximum  distance  no  contribution  is  made  to  the  Green's  function  integral.  It  is 
not  expected  that  this  will  be  possible  for  objects  less  than  a  few  wavelengths  in 
maximum  dimension. 

•  Modify  the  mesh  generation  program  to  allow  for  conformity  to  multiple  interfaces 
(muiti-layered  objects  without  the  granular  noise  problem). 

•  Modify  all  programs  for  Intel  80386  based  Fortran  compiler  use,  such  as  NDP 
FORTRAN  by  Microway.  This  will  remove  the  640  KB  memory'  limitation  im¬ 
posed  by  the  IBM  DOS  and  allow's  for  the  solution  to  larger  scattering  problems. 
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APPENDIX  A.  INPUT  DATA  FILE  EXAMPLE 


CIRCLE 

R 

NI 

ND 

NU 

NM 

0.019 
0.  03 
1.0 
0 
0 

999.  9 
36 
72 
1 
5 

35 
4.5 
4.  0 
0.  3 
1.  0 
0.  0 
1.  0 
0.  0 
P 

1.0 

3.  OOE+08 


-  PLOT  TITLE 

-  RECTANGULAR  INPUT  DATA 

-  NO  INTERMEDIATE  DATA  RECORDING 

-  NO  DISSPLA  PROGRAM  GENERATION 

-  NON  UNIFORM  MATERIAL  (MATERIAL  INTERFACE  PRESENT) 

-  NO  MESH  GENERATION  ONLY 

-  REQUESTED  MESH  RESOLUTION  (IN  WAVELENGTHS) 

-  CONTOUR  DISTANCE  (IN  WAVELENGTHS) 

-  MULTIPLICATIVE  GFI  SCALE  FACTOR 

-  DISTANCE  <1.0  BIAS  TERM 

-  DISTANCE  >  1.  0  BIAS  TERM 

-  MAX  DIST  BEYOND  WHICH  THERE  IS  NO  CONTRIBUTION  TO  GFI 

-  NUMBER  OF  INPUT  DATA  POINTS 

-  NUMBER  OF  POINTS  FOR  SIGMA  CALCULATION  (IN  THE  CIRCLE) 

-  ELEMENT  GENERATION  METHOD 

-  START  NODE 

-  STOP  NODE 

-  X  AXIS  OFFSET 

-  Y  AXIS  OFFSET 

-  DESIRED  DIMENSION  (ORIGIN  TO  FIRST  POINT  .WAVELENGTHS) 

-  REAL  PART  OF  ALPHA 

-  IMAGINARY  PART  OF  ALPHA 

-  REAL  PART  OF  BETA 

-  IMAGINARY  PART  OF  BETA 

-  PLANEWAVE  GENERATOR  ENABLED 

-  PLANEWAVE  AMPLITUDE 

-  INPUT  FREQUENCY  OF  THE  PLANEWAVE,  FO 


0 

. 00000000 

. 50000000 

10 

. 08682409 

. 49240390 

20 

. 17101010 

. 46984630 

30 

. 25000000 

.43301270 

40 

. 32139380 

. 38302220 

50 

. 38302220 

. 32139380 

60 

. 43301270 

. 25000000 

70 

. 46984630 

. 17101010 

80 

. 49240390 

. 08682407 

90 

. 50000000 

-. 00000002 

100 

.49240390 

-.08682411 

110 

. 46984630 

-. 17101010 

120 

.43301270 

-.25000000 

130 

. 38302220 

-. 32139380 

140 

. 32139380 

-. 38302220 

150 

. 25000000 

-.43301270 

160 

. 17101000 

-. 46984630 

170 

. 08682404 

-.49240390 

ISO 

-. 00000004 

-. 50000000 

190 

-.08682413 

-.49240390 

200 

17101010 

-.46984630 

210 

-. 25000000 

-.43301270 

220 

-. 32139380 

-. 38302220 

230 

-. 38302220 

-. 32139380 

-  ANGLE,  X,  Y  COORD. 
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240 

-.43301270 

250 

-.46984630 

260 

-.49240390 

270 

-. 50000000 

280 

-. 49240390 

290 

-.46984630 

300 

-.43301270 

310 

-. 38302220 

320 

-. 32139380 

330 

-. 24999990 

340 

-. 17101000 

350 

-.08682401 

25000000 
17101000 
08682403 
.  00000007 
.  08682416 
.  17101010 
.  25000010 
.  32139390 
.  38302230 
.43301280 
.46984630 
.49240390 
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APPENDIX  B.  READ  PROGRAM 


REAL  A(0:2000),  B(0:2000),  MRES,  DIST,  XORIGIN,  YORIGIN 
REAL  C,  D,  E,  F,  DPER,  DD,  FREQ,  EO,  MAXD 

INTEGER  I,  J,  METHOD,  STARTND,  STOPND,  NRES,  MODE,  LB IAS,  GBIAS 
CHARACTER* 1  CHAR, CHAR 1,CHAR2,CHAR3, CHAR4 , CHARS 
CHARACTER* 13  NAME 

OPENCUNIT  =  1,  FILE  =  'D:  FINALDWG.DAT’) 

OPENCUNIT  =  2,  FILE  =  'C:  MSFORT  INPUT.DAT') 

J  =  0 

WRITEC*,*)  'ENTER  THE  PLOT  NAME  OR  LABEL  (MAX  OF  13  CHARACTERS)’ 
READ(*,1005)  NAME 

WRITEC*,*)  'ENTER  THE  COORDINATE  SYSTEM  IN  USE.  (R  OR  P)' 
READ(*,1000)  CHAR 

WRITEC*,*)  'DO  YOU  WANT  INTERMEDIATE  VALUES  ?  (I)' 

READ(*,1000)  CHARI 

WRITEC*,*)  'DO  YOU  WANT  A  DISSPLA  PROGRAM  GENERATED  ?  (D)' 
READ(*,1000)  CHAR2 

WRITEC*,*)  'UNIFORM  SLAB  (NO  MATERIAL  TO  VACUUM  INTERFACE  ?  (U)' 
READ(*,1000)  CHAR3 

WRITEC*,*)  'MESH  GENERATION  ONLY  ?  (M)’ 

READ(*,1000)  CHAR4 

WRITEC*,*)  'ENTER  THE  MESH  RESOLUTION.  (0.4,  ETC...)’ 

RE ADC*,*)  MRES 

WRITEC*,*)  'ENTER  THE  CONTOUR  DISTANCE.  (0.3,  ETC...)’ 

READC*,*)  DIST 

WRITEC*,*)  'ENTER  THE  GFI  SCALE  FACTOR. Cl. 0,  ETC...)' 

READC*,*)  DPER 

WRITEC*,*)  'ENTER  THE  GFi  C  <  1  BIAS  TERM  C0,1,  ETC...)' 

READC*,*)  LB IAS 

WRITEC*,*)  'ENTER  THE  GFI  C  >  1  BIAS  TERM  C0,1,  ETC...)' 

READC*,*)  GBIAS 

WRITEC*,*)  'ENTER  THE  GFI  MAX  DISTANCE  (999.0,  ETC...)' 

READC*,*)  MAXD 

WRITEC*,*)  'ENTER  THE  NUMBER  OF  POINTS  FOR  SIGMA  CALC. C36,ETC.  .  . ) ' 
READC*,*)  NRES 

WRITEC*,*)  'ENTER  THE  DRAWING  METHOD.  Cl  -  6)' 

READC*,*)  METHOD 

WRITEC*,*)  'ENTER  THE  STARTING  NODE  NUMBER.  C MUST  BE  >  0)' 

READC*,*)  STARTND 

WRITEC*,*)  'ENTER  THE  BISECTION  STOPPING  NODE  NUMBER.  ' 

READC*,*)  STOPND 

WRITEC*,*)  'DESIRED  DISTANCE  FROM  ORIGIN  TO  POINT  1  CIN  WLs)  ' 
READC*,*)  DD 

WRITEC*,*)  'ENTER  THE  X  AXIS  ORIGIN. ’ 

READC*,*)  XORIGIN 

WRITEC*,*)  'ENTER  THE  Y  AXIS  ORIGIN. ' 

READC*,*)  YORIGIN 

WRITEC*,*)  'ENTER  THE  REAL  COMPONENT  OF  ALPHA. ' 

READC*,*)  C 
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WRITE(*,*)  'ENTER  THE  IMAGINARY  COMPONENT  OF  ALPHA.  ’ 
READ(*,*)  D 

WRITE(*,*)  'ENTER  THE  REAL  COMPONENT  OF  BETA.  ' 
READ(*,*)  E 

WRITE(*,*)  'ENTER  THE  IMAGINARY  COMPONENT  OF  BETA.  ' 
READ(*,*)  F 

WRITE(*,*)  'plane  wave  (P)  OR  CYLINDRICAL  MODE  (C)  ' 
READ(*,1000)  CHARS 

WRITE(*,*)  'ENTER  THE  WAVE  FREQUENCY  (IN  HERTZ)' 
READ(*,*)  FREQ 

WRITE(*,*)  'ENTER  THE  WAVE  AMPLITUDE  ’ 

READ(*,*)  EO 

IF((CHAR5.EQ.  'C').OR.  (CHARS.  EQ.  'c'))  THEN 
WRITE(*,*)  ’enter  the  MODE  NUMBER  ' 

READ(*,*)  MODE 

ENDIF 

DO  10,  1=0,  1999 

READ(1,*,ERR  =  20)  A(J),  B(J) 

IF(A(J).  GT.  1000)  THEN 
GOTO  S 

ELSEIF(B(J).GT.  1000)  THEN 
GOTO  S 

ELSE 

J  =  J  +  1 

ENDIF 

S  CONTINUE 

10  CONTINUE 
20  J  =  J  -  1 

WRITE(2,100S)  NAME 
WRITE(2,1000)  CHAR 
WRITE(2,1000)  CHARI 
WRITE(2,1000)  CHAR2 
WRITE(2,1000)  CHARS 
WRITE(2,1000)  CHAR4 
WRITE(2,1010)  MRES 
WRITE(2,1010)  DIST 
WRITE(2,1010)  DPER 
WRITE(2,1020)  LBIAS 
WRITE(2,1020)  GBIAS 
WRITE(2,1040)  MAXD 
WRITE(2,1020)  J 
WRITE(2,1020)  NRES 
WRITE(2,1020)  METHOD 
WRITE(2,1020)  STARTND 
WRITE(2,1020)  STOPND 
WRITE(2,1010)  XORIGIN 
WRITE(2,1010)  YORIGIN 
WRITE(2,1010)  DD 
WRITE(2,1010)  C 
WRITE(2,1010)  D 
WRITE(2,1010)  E 
WRITE(2,1010)  F 
WRITE(2,1000)  CHARS 
WRITE(2,1010)  EO 

IF( (CHARS.  EQ.  ' C ).  OR.  (CHARS.  EQ.  'c'))  THEN 
WRITE(2,1020)  MODE 
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ENDIF 

WRITE(2,*)  FREQ 
DO  30,  1=1,  J 

WRITE(2,1030)  I,  A(I),  B(I) 
30  CONTINUE 

WRITE(2,*)  'END  OF  FILE' 

C 

1000  FORMAT(Al) 

1005  F0RMAT(A13) 

1010  F0RMAT(F8. 6) 

1020  F0RMAT(I3) 

1030  F0RMAT(1X,I3,7X,F12.  8,5X,F12.  8) 
1040  F0RMAT(F8,  3) 

C 

CLOSE(l) 

CL0SE(2) 

STOP 

END 

C 
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APPENDIX  C.  MESH  GENERATION/FINITE  ELEMENT  PROGRAM 


FINITE  ELEMENT  MESH  GENERATION  PROGRAM 
WRITTEN  BY  LT.  T.  B.  WELCH 

w/  PROGRAMMING  IDEAS  FROM  PROF.  M. A.  MORGAN 

COMPLETELY  FILE  (INPUT.DAT)  DRIVEN 

INPUT  PARAMETERS: 

CHARACTER  (POLAR  OR  RECTANGULAR  INPUT  DATA  -  FLAG) 
CHARACTER  (INTERMEDIATE  MATRICES  WRITE  -  FLAG) 

CHARACTER  (DISSPLA  PROGRAM  GENERATION  -  FLAG) 

CHARACTER  (UNIFORM  MATERIAL  -  FLAG) 

CHARACTER  (MESH  GENERATION  ONLY  -  FLAG) 

MESH  RESOLUTION 
CONTOUR  OFFSET 

DESIRED  SCALE  FACTOR  FOR  GREEN’S  FUNCTION  INTEGRAL 
BIAS  (SHIFT)  FOR  GFI  STEP  WHEN  <1.0 
BIAS  (SHIFT)  FOR  GFI  STEP  WHEN  >1.0 
MAXIMUM  DISTANCE  BEYOND  WHICH  THERE  WILL  BE  NO 

CONTRIBUTION  TO  THE  GREEN’S  FUNCTION  INTEGRAL 
NUMBER  OF  POINTS 

NUMBER  OF  POINTS  FOR  CROSS  SECTION,  (IN  THE  CIRCLE) 

MESH  GENERATION  TECHNIQUE 
START  NODE 
STOP  NODE 
X  ORIGIN 
Y  ORIGIN 

DESIRED  DISTANCE  (FROM  ORIGIN  TO  FIRST  POINT,  WAVELENGTH) 
ALPHA  (INPUT  AS  A  AND  B  THEN  CONVERTED  TO  COMPLEX  ALPHA) 
BETA  (INPUT  AS  A  AND  B  THEN  CONVERTED  TO  COMPLEX  BETA) 
CHARACTER  (PLANEWAVE  GENERATOR  -  FLAG) 

AMPLITUDE 

FREQUENCY  (IN  HERTZ) 

--  OR  -- 

AMPLITUDE 

CYLINDRICAL  MODE  NUMBER  (GENERATOR  ALWAYS  DOES  N=l) 
FREQUENCY  (IN  HERTZ) 

--  OR  -- 

BOUNDARY  CONDITIONS 

DATA  POINT  PAIRS  (X,  Y  OR  RADIUS,  THETA) 

END  OF  FILE  MARKER 


COMPLEX  ALPHA, BETA, BCOND( 100) ,LINE(50) ,ANS( 100) 

COMPLEX  A(50,50),B(50,50),C(50,50),P(50,100),SURBC(100) 

COMPLEX  U(100,100),PSI(50,100),SVEC(50,100) 

REAL  MRES,MESH(0:  200,5) ,NDP(200,2) ,OFFSET,DPER,MRESW 
REAL  MINAREA .MAXAREA , AREA ,E0 ,XORIGIN, YORIGIN ,K0 ,MAXD 
INTEGER  NPNTS , PERND , N0DES( 200 ) , MOR( 200 ) , B IND , NBMX , UNK , LB I AS , GB I AS 
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c 


c 


c 


c 


c 


INTEGER  LND( 0:200,3), NDL( 200 , 4) , NCT( 200 ) , MAXEL , IMX , MODE , NRES 
INTEGER  METHOD , STARTND , STOPND , MINROW , MINED , MAXROW , NABC (100,3) 
CHARACTER*!  CHAR,  CHARI,  CHAR2,  CHAR3,  CHAR4,  CHARS 
EQUIVALENCE  (PS I,  P) 


COMMON/BLKl/MESH 
C0MM0N/BLK2/PERND , BIND 
C0MM0N/BLK3/N0DES 
C0MM0N/BLX4/LND ,NDL ,NCT 
C0MM0N/BLK5/NDP 
C0MM0H/BLK6/M0R 

C0MM0N/BLK7 /MINAREA . MINROW , MINED ,MAXAREA , MAXROW , MAXEL , AREA 

C0MM0N/BLK8/A,B,C,P 

COMMON/ BLK9/CHAR2,  CHAR3,  CHAR4 


OPEN  (UNIT  = 
OPEN  (UNIT  = 
OPEN  (UNIT  = 
OPEN  (UNIT  = 
OPEN  (UNIT  = 


1,  FILE  = 

2,  FILE  = 

3,  FILE  = 

4,  FILE  = 
7,  FILE  = 


CACCESS=' SEQUENTIAL' ,FORM= 
OPEN  (UNIT  =  8,  FILE  =  ’ 
CACCESS=’ SEQUENTIAL' ,FORM^ 


OPEN 

(UNIT 

= 

9, 

FILE 

— 

OPEN 

(UNIT 

10, 

FILE 

OPEN 

(UNIT 

11, 

FILE 

OPEN 

(UNIT 

12, 

FILE 

= 

OPEN 

(UNIT 

S 

13, 

FILE 

= 

OPEN 

(UNIT 

ss 

19, 

FILE 

= 

OPEN 

(UNIT 

= 

20, 

FILE 

S 

OPEN 

(UNIT 

= 

30, 

FILE 

S= 

OPEN 

(UNIT 

s 

40, 

FILE 

= 

MINAREA  =  999999, 

,  9 

MAXAREA  =  • 

•999999. 9 

MATRIX.  DAT’ , STATUS  =  'UNKNOWN') 
PLT.  DAT' , STATUS  =  'UNKNOWN') 
INPUT.DAT' , STATUS  =  'OLD') 

TEXT.  LBL' , STATUS  =  'UNKNOWN' ) 
RS.DAT' , STATUS  =  'UNKNOWN' , 
UNFORMATTED') 

PSI.  DAT’ , STATUS  =  'UNKNOWN' , 
FORMATTED'  ) 

ABCP.DAT' .STATUS  =  'UNKNOWN') 
BCOND.dat',  status  =  'UNKNOWN') 
FMATRIX.  DAT' ,STATUS=' UNKNOWN' ) 
NABC.  DAT’ ,STATUS=' UNKNOWN' ) 
OUTPUT.  DAT' ,STATUS=’ UNKNOWN' ) 
ABRMAT.  DAT' , STATUS= ’ UNKNOWN ' ) 
DISPLA.  DAT' ,STATUS=' UNKNOWN' ) 

U.  DAT' ,STATUS=' UNKNOWN' ) 

F3.  DAT’ ,STATUS=' UNKNOWN’ ) 


CALL  IO( NPNTS , MRES , METHOD , STARTND , STOPND , OFFSET , ALPHA , BETA , BCOND , 
CE 0 , CHARI , MODE , XORI GIN , YORI GIN , CHARS , DPER , KO , NRES , MRESW , LB IAS , GB I AS 
C,MAXD) 

CALL  ROTATE ( NPNTS , METHOD , STARTND , STOPND ) 

CALL  BOUND( MRES, NPNTS, METHOD, STOPND) 

CALL  NORMAL 

CALL  NODSET( METHOD, NABC, NBMX,UNK) 

PAUSE  'PLEASE  PRESS  ENTER  TO  CONTINUE,  OR  CTRL  BREAK  TO  ABORT!  ' 
CALL  LOADER( BCOND , OFFSET , ALPHA , BETA , NABC , IMX , NBMX , E 0 , SURBC , CHARI , 
CLINE ,MODE , XORIGIN , YORIGIN , SVEC , CHARS ) 

CALL  SWEEP( IMX , NABC , SURBC , LINE , CHAR 1 , U , BCOND , PS I , ANS , CHARS ) 

CALL  SAVE( BCOND , ANS , U , OFFSET, PERND , CHARS , DPER , KO , XORIGIN , YORIGIN , 
CNRES , MRESW , LB IAS , GB IAS , MAXD ) 

CLOSE(l) 

CLOSE(2) 

CLOSE(3) 

CL0SE(4) 

CLOSE(7) 

CLOSE(8) 
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CLOSEO) 

CLOSE(IO) 

CLOSE(ll) 

CLOSE(12) 

CLOSE(13) 

CLOSE(19) 

CL0SE(20) 

CL0SE(30) 

CLOSE(40) 

STOP 

END 

SUBROUTINE  IO( NPNTS , MRES , METHOD , STARTND , STOPND , DIST , ALPHA , BETA , 
CBCOND , EO , CHARI , MODE , XORIGIN , YORIGIN , CHARS , DPER , KO , NRES , MRESW , 
CLBIAS,GBIAS,MAXD) 

THIS  SUBROUTINE  READS  IN  THE  INPUT  PARAMETERS 
AND  SURFACE  DATA  POINTS.  THESE  POINTS  CAN  BE 
IN  EITHER  POLAR  OR  RECTANGULAR  FORM. 

COMPLEX  ALPHA, BETA, BCOND( 100) 

REAL  A , B , PI , XORIGIN , YORIGIN , EO ,KO , FO , SFAC , DD 

REAL  MESH(0:  200,5) , MRES, THETA, RADIUS, DIST, C, LAMBDA, MRESW 

REAL  DISTW,DPER,MAXD 

INTEGER  I , I I , J , K , NPNTS , RES , METHOD , STARTND , STOPND , MODE , NRES , LB I AS 
INTEGER  GBIAS 

CHARACTER* 1  CHAR , CHAR 1 , CHAR2 , CHAR3 , CHAR4 , CHARS 
CHARACTER* 12  NAME 

COMMON/BLKl/MESH 
COMMON/BLK9/CHAR2,  CHARS,  CHAR4 

C  =  2. 997925E+08 
PI  =  4. 0*ATAN(1. 0) 

READ(3,1070)  NAME 
READ(3,1000)  CHAR 
READ(3,1000)  CHAR2 
READ(3,1000)  CHARS 
READ(3,1000)  CHAR4 
READ(3,1000)  CHARS 
READ(3,*)  MRESW 
READ(3,*)  DISTW 
READ(3,’)  DPER 
READ(3,*)  LBIAS 
READ(3,*)  GBIAS 
READ(3,*)  MAXD 
READ(3,*)  RES 
READ(3,*)  NRES 
READ(3,*)  METHOD 
READ(3,*)  STARTND 
READ(3,*)  STOPND 
READ(3,*)  XORIGIN 
READ(3,*)  YORIGIN 
READ(3,*)  DD 
READ(3,*)  A 
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READ(3,*)  B 
ALPHA  =  CMPLX(A,B) 

READ(3,*)  A 
READ(3,*)  B 
BETA  =  CMPLX(A,B) 

READ(3,1000)  CHARI 
C 

IFCMETHOD.EQ.  1)  THEN 

WRITE(6,*)  'METHOD  1  SELECTED  -  OBJECT  BISECTION  ' 
ELSEIFCMETHOD.EQ.  2)  THEN 

WRITE(6,*)  'METHOD  2  SELECTED  -  CONNECTED  NODES  ’ 
ELSEIFCMETHOD.EQ.  3)  THEN 

WRITE(6,*)  'METHOD  3  SELECTED  -  SELECTED  STOP  NODE  ' 
ELSEIFCMETHOD.EQ.  4)  THEN 

WRITEC6,*)  'METHOD  4  SELECTED  -  CONNECTED/ SELECTED  STOP  NODE 
ELSEIFCMETHOD.EQ.  5)  THEN 

WRITEC6,*)  'METHOD  5  SELECTED  -  EQUALLY  SPACED  CONNECTD  NODE 

ELSE 

WRITEC6,*)  'METHOD  6  SELECTED  -  SELECTED  START  AND  STOP  NODE 
ENDIF 
C 

IFCCCHAR2.  EQ.  'l').OR.  CCHAR2.  EQ.  'i'))  THEN 

WRITEC6,*)  'INTERMEDIATE  MATRIX  FILE  GENERATION  -  ENABLED' 

ELSE 

WRITEC6,*)  'INTERMEDIATE  MATRIX  FILE  GENERATION  -  DISABLED 

ENDIF 

C 

IFCCCHAR3.EQ.  'D').OR.  CCHAR3.EQ.  'd'))  THEN 

WRITEC6,*)  'DISSPLA  FORTRAN  PROGRAM  GENERATION  -  ENABLED' 

ELSE 

WRITEC6,*)  'DISSPLA  FORTRAN  PROGRAM  GENERATION  -  DISABLED 

ENDIF 

C 

IFCCCHAR4.EQ. 'U').OR. CCHAR4.EQ. 'u'))  THEN 

WRITEC6,*)  'UNIFORM  MATERIAL  SPECIFIED  CNO  INTERFACE)  ' 

ELSE 

WRITEC6,*)  'MATERIAL  SPECIFIED  WITH  A  VACUUM  AROUND  OBJECT' 

ENDIF 

C 

IFCCCHAR5.EQ.  'M').0R.  CCHAR5.EQ. 'm'))  THEN 

WRITEC6,*)  'MESH  GENERATION  «<  ONLY  »>  -  ENABLED' 

ELSE 

WRITEC6,*)  'MESH  GENERATION  AND  FE  PROGRAM  -  ENABLED' 

ENDIF 

C 

IFCCCHARl.EQ.  'P').OR.  CCHARl.EQ.  'p'))  THEN 
READC3,*)  EO 
READC3,*)  FO 
LAMBDA  =  C/FO 
KO  =  2*PI /LAMBDA 

WRITEC6,*)  'PLANEWAVE  BOUNDARY  VALUE  GENERATION  -  ENABLED' 
WRITEC6,*)  'AMPLITUDECEO)  =  EO,',  WAVENUMBERCKO)  =  ' ,KO 

WRITEC 6,*) 'WAVELENGTH  =' .LAMBDA, ',  FREQUENCYCFO)  =  ' ,F0 
ELSEIFC CCHARl.EQ.  'C').OR.  CCHARl.EQ.  'c'))  THEN 
READC3,*)  EO 
READC3,*)  MODE 
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READ(3,*)  FO 
LAMBDA  =  C/FO 
KO  =  2*PI*LAMBDA 

WRITE(6,*)  'CLYINRICAL  BOUNDARY  VALUE  GENERATION  -  ENABLED' 
WRITE(6,*)  'AMPLITUDE(EO)  =  EO,',  MODE  NUMBER  =  MODE 

WRITE( 6,*) 'WAVELENGTH  =' .LAMBDA, ',  FREQUENCY(FO)  =  ' ,FO 

ELSE 

WRITE(6,*)  'ALL  BOUNDARY  VALUE  GENERATION  METHODS  -  DISABLED' 
DO  5,  K  =  1,  RES 

READ(3,*)  BCOND(K) 

5  CONTINUE 

ENDIF 

POLAR  COORDINATE  INPUT  ROUTINE 

IF((CHAR.EQ.  'P').OR.  (CHAR.EQ.  'p'))  THEN 
READ(3,1020)  THETA,  RADIUS 
SFAC  =  2*PI*DD/RADIUS 

MESH(0,4)  =  (RADIUS*SIN(THETA*PI/180.  0))*SFAC  +  XORIGIN 
MESH(0,5)  =  (RADIUS*COS(THETA*PI/180.  0))*SFAC  +  YORIGIN 
DO  10,  J  =  1,  RES-1 

READ(3,1020)  THETA,  RADIUS 

MESH(J,4)  =  (RADIUS*SIN(THETA*PI/180.  0))*SFAC  +  XORIGIN 
MESH(J,5)  =  (RADIUS*COS(THETA*PI/180. 0))*SFAC  +  YORIGIN 
CONTINUE 

WRITE(6,*)  'POLAR  COORDINATE  INPUT  SELECTED  ' 

RECTANGULAR  COORDINATE  INPUT  ROUTINE 

ELSEIF(( CHAR.EQ.  ' R ' ) .  OR.  ( CHAR. EQ.  'r'))  THEN 
READ(3,1010)  MESH(0,4),  MESH(0,5) 

SFAC  =  2*PI*DD/((MESH(0,4)**2+MESH(0,5)**2)**0. 5) 
MESH(0,4)  =  MESH(0,4)*SFAC  +  XORIGIN 

MESH(0,5)  =  MESH(0,5)*SFAC  +  YORIGIN 

DO  20,  J  =  1,  RES‘l 

READ(3,1010)  MESH(J,4),  MESH(J,5) 

MESH(J,4)  =  MESH(J,4)*SFAC  +  XORIGIN 

MESH(J,5)  =  MESH(J,5)*SFAC  +  YORIGIN 

20  CONTINUE 

WRITE(6,*)  'RECTANGULAR  COORDINATE  INPUT  SELECTED  ' 

ELSE 

WRITE(6,*)  'INPUT  DATA  FILE  COORDINATE  SPECIFICATION  ERROR' 

ENDIF 

WRITE(*,1100)  DD,  SFAC 
MESH(RES,4)  =  MESH(0,4) 

MESH(RES,5)  =  MESH(0,5) 

NPNTS  =  J 

WRITE(6,*)  'NODE  AT  0  DEGREES  (X/Y  COORDINATES)  =  ',  MESH(0,4), 
CMESH(0,5) 

WRITE(6,*)  X,  Y  OFFSETS  =  ',  XORIGIN,  YORIGIN 
WRITE(6,1090)  ALPHA,  BETA 
MRES  =  MRESW*2*PI 

WRITE(6,*)  'MESH  RESOLUTION  =  '.MRESW,'  WAVELENGTHS' 

DIST  =  DISTW*2*PI 

WRITE(6,*)  'CONTOUR  DISTANCE  =  ',DISTW,'  WAVELENGTHS' 

WRITE(6,*)  'NUMBER  OF  DATA  POINTS  =  ' ,RES 
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IF((METHOD,EQ.  3).  OR.  (METHOD.  EQ.  6))  THEN 

WRITE(6,*)  'START  NODE  =  ' .STARTND 

WRITE(6,*)  'STOP  NODE  =  ' .STOPND 


ENDIF 

WRITE(4,*) 

'17' 

WRITE(4,*) 
C  55 

'10 

1' 

1 

0 

WRITE(4,*) 
C  80 

'2 

1' 

1 

0 

WRITE(4,*) 
C  90 

'2 

2' 

1 

0 

WRITE(4,*) 
C  100 

'2 

1' 

1 

0 

WRITE(4,*) 
C  110 

'2 

2' 

1 

0 

WRITE(4,*) 
C  120 

'2 

1' 

1 

0 

WRITE(4,*) 
C  130 

'2 

2' 

1 

0 

WRITE(4,*) 

C  140 

'2 

1' 

1 

0 

WRITE(4,*) 
C  150 

'2 

2' 

1 

0 

WRITE(4,*) 

C  160 

'2 

1' 

1 

0 

WRITE(4,*) 
C  170 

'2 

2' 

1 

0 

WRITE(4,*) 
C  180 

'2 

1' 

1 

0 

WRITE(4,*) 
C  190 

'2 

2' 

1 

0 

WRITE(4,*) 
C  200 

'2 

1' 

1 

0 

WRITE(4,*) 
C  210 

'2 

2' 

1 

0 

WRITE(4,*) 
C  220 

'2 

1' 

1 

0 

WRITE(4,*) 
C  230 

'2 

2' 

1 

0 

VRITE(4,1070)  NAME 
WRITE(4,*)  'L' 

VfRITE(4,*)  'Mesh  Resolution' 
WRITE(4,*)  'L' 

WRITE(4,1030)  MRESW 
WRITE.(4,*)  'L' 

WRITE(4,*)  'Contour  Distance' 
WRITE(4,*)  'L' 

WRITE(4,1030)  DISTW 
WRITE(4,*)  'L' 

WRITE(4,*)  'Number  of  Points' 
WRITE(4,*)  'L' 

WRITE(4,1040)  RES 
WRITE(4,*)  'L' 

WRITE (4,*)  'Method' 

WRITE(4,*)  'L' 


130 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 
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WRITE(4,1040)  Method 
WRITE(4,*)  ’L’ 
WRITE(4,*)  'X  Origin' 
WRITE(4,*)  'L' 
WRITE(4,1050)  XORIGIN 
WRITE(4,*)  'L' 
¥RITE(4,*)  'Y  Origin' 
WRITE(4,*)  'L' 
WRITE(4,1050)  YORIGIN 
WRITE(4,*)  'L' 
WRITE(4,*)  'Alpha' 
WRITE(4,*)  'L' 
WRITE(4,1060)  ALPHA 
WRITE(4,*)  'L' 
WRITE(4,*)  'Beta' 
WRITE(4,*)  'L' 
WRITE(4,1060)  BETA 
WRITE(4,*)  'L' 
WRITE(4,*)  '9' 


30 


C 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1090 

1100 


DO  30,  II  =  0,  NPNTS-1 

WRITE(2,*)  MESH(II,4),  MESH(II.5) 
CONTINUE 

WRITE(2,*)  MESH(0,4),  MESH(0,5) 

WRITE(2,*)  '999992,  999991  ' 

WRITE(2,*)  '999990,  999990  ' 

FORMAT(Al) 

FORMAT(10X,F12.  8,4X,F12.  8) 
FORMAT(3X,I3,5X,F12,  8) 

F0RMAT(F8.  6) 

F0RMAT(I4) 

FORMATC 1X,F8. 6) 

F0RMAT(1X,F8.6,'  -  J' ,F8. 6) 

F0RMAT(A12) 

FORMATC/) 

FORMATC IX, 'ALPHA  =  ',F8.4,2X,'+  J  ',Fe.4.' 
C  '  ,F6.4) 

FORMATC IX. 'DESIRED  DISTANCE  =  ',F8.5,'. 

RETURN 

END 


BETA  =  'F8.  4,2X,'+  J 
SCALE  FACTOR  =  ',F8. 5) 


SUBROUTINE  ROTATE C  NPNTS , METHOD , STARTND , STOPND ) 

THIS  SUBROUTINE  ROTATES  THE  SURFACE  POINTS  TO  ALLOW  FOR 
A  REARRANGING  OF  THE  START  NODE  CFIRST  DATA  POINT). 

REAL  MESHCO: 200,5) 

INTEGER  I,  NPNTS,  METHOD,  STARTND,  STOPND 
COMMON/BLKl/MESH 


IFCMETHOD.  EQ.  6)  THEN 

DO  10,  I  =  0,  NPNTS-1 

MESHCI.l)  =  MESHCI,4) 
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MESH(I,2)  =  MESH(I,5) 

10  CONTINUE 

C 

DO  20,  I  =  STARTND-1,  NPNTS-1 

MESH(I-(STARTND-1),4)  =  MESH(I,1) 
MESH(I-( STARTND-1), 5)  =  MESH(I,2) 

D  CONTINUE 

DO  30,  1=1,  STARTND 

MESH(I+NPNTS-STARTND,4)  =  MESH(I-1,1) 
MESHd+NPNTS -STARTND,  5)  =  MESH(  1-1,2) 
0  CONTINUE 

ENDIF 

STOPND  =  STOPND  -  STARTND  +  1 
IF(ST0PND.  GT.NPNTS)  THEN 

STOPND  ■=  STOPND  -  NPNTS 
ELSEIFCSTOPND.LT.  0)  THEN 

STOPND  =  NPNTS  +  STOPND  +1 
ELSEIFCSTOPND.EQ.  0)  THEN 
STOPND  =  1 

ELSE 

STOPND  =  STOPND 

ENDIF 

RETURN 

END 


SUBROUTINE  BOUND( MRES , NPNTS , METHOD , STOPND) 

THIS  SUBROUTINE  REORDERS  THE  THE  SURFACE  POINTS 
BASED  ON  THE  DESIRED  INPUT  MESH  RESOLU  INN.  THE 
OBJECT  IS  ALSO  BISECTED. 

REAL  PERIM,  DIST,  MESH(0: 200 ,5) ,  TEMP,  TEMPI,  MRES,  MRESN,  MRESLW 
REAL  MRESNL,  MRESNR,  TEMPR,  TEMPL,  D.TSTB,  DZ,  PI,  MRESW,  MRESRW 
INTEGER  I,  PERND,  BIND,  NPNTS , STARTND ,  STOPND,  METHOD,  NEWND,  J 
INTEGER  M,  L 
COMMON/BLKl/MESH 
COMMON/BLK2/PERND , B IND 
C 

PI  =  4. 0*ATAN(1.  0) 

PERIM  =  0.  0 
C 

DO  10,  1=0,  NPNTS-1 

DIST  =  SQRT((MESH(I+1,4)-MESH(I,4))**2+(MESH(I+1,5)  - 
C  MESH(I,5))**2) 

PERIM  =  PERIM  +  DIST 
MESH( 1+1,3)  =  PERIM 
10  CONTINUE 

WRITE(6,*)  'PERIMETER  LENGTH  =  ',  PERIM 

PERND  =  NINT(PERIM/MRES  -  AMOD(PERIM/MRES,2.  0)) 

BIND  =  (PERND  -  2)/2 

WRITE(6,*)  'PERIMETER  NODE  =  ',  PERND 

MRESW  =  MRES/(2*PI) 

WRITE(6,*)  'YOUR  REQUESTED  MESH  RESOLUTION  OF  ', MRESW 
IF( (METHOD.  EQ.  3).  OR.  (METHOD.  EQ.  4).  OR.  (METHOD.  EQ.  5).  OR. 


84 


o  o  o 


CCMETHOD.EQ.  6))  THEN 

MRESNL  =  (PERIM  -MESH(STOPND-l,3))/FLOAT(PERND/2) 

MRESLW  =  MRESNL/(2*PI) 

MRESNR  =  MESH(STOPND-l,3)/FLOAT(PERND/2) 

MRESRW  =  MRESNR/ (2*PI) 

WRITE(6,*)  .  HAS  BEEN  MODIFIED  TO  .  .  .  ' 

WRITE(6,*)  'LEFT  SIDE  OF  SEGMENT  .  .  .  ’.MRESLW 
WRITE(6,*)  'RIGHT  SIDE  OF  SEGMENT  .  .  .  MRESRW 

ELSE 

MRESN  =  PERIM/FLOAT(PERND) 

MRESW  =  MRESN/ (2*PI) 

WRITE(6,*)  .  HAS  BEEN  MODIFIED  TO  .  .  .  '.MRESW 

ENDIF 

PERIMETER  NODE  INITIALIZATION  (X.Y  COORD) 

IF( (METHOD.  EQ.  1). OR. (METHOD. EQ. 2))  THEN 
DO  30.  1=0.  PERND-1 
TEMP  =  MRESN*I 
J  =  1 

20  IF(TEMP.  GT.MESH(J.3))  THEN 

J  =  J  +  1 
GOTO  20 

ENDIF 

TEMPI  =  (TEMP  -  MESH(J-1.3))/(SQRT((MESH(J,4)-MESH(J-1.4)) 

C  **2  +  (MESH(J.5)  -  MESH(J-1.5))**2)) 

MESH(I+1,1)  =  MESH(J-1.4)  +  TEMP1*(MESH( J.4)-MESH( J-1.4)) 
MESH(I+1.2)  =  MESH(J-1.5)  +  TEMP1*(MESH( J.5)-MESH( J-l.S)) 

30  CONTINUE 

C 

ELSE 

DO  60,  1=0.  PERND-1 
TEMPR  =  MRESNR*! 

TEMPL  =  PERIM  -  MRESNL*( PERND  -  I) 

IF(TEMPR.  LE.  MESH(STOPND-l .3) )  THEN 
J  =  1 

40  IF(TEMPR.  GT.MESH(J,3))  THEN 

J  =  J  +  1 
GOTO  40 

ENDIF 

TEMPI  =  (TEMPR  -  MESH(J-1.3))/(SQRT((MESH(J.4)-MESH( J-1.4)) 
C  **2  +  (MESH(J.5)  -  MESH(J-1.5))**2)) 

MESH(I+1.1)  =  MESH( J-1.4)  +  TEMP1*(MESH(J,4)-MESH( J-1.4)) 
MESH(I+1,2)  =  MESH(J-1.5)  +  TEMP1*(MESH(J.5)-MESH(J-1.5)) 
ELSE 

J  =  1 

50  IF(TEMPL. GT.MESH(J.3))  THEN 

J  =  J  +  1 
GOTO  50 

ENDIF 

TEMPI  =  (TEMPL  -  MESH(J-1.3))/(SQRT((MESH(J.4)-MESH( J-1.4)) 
C  **2  +  (MESH(J,5)  -  MESH(J-1,5))**2)) 

MESH(I+1.1)  =  MEoIK  J-1.4)  +  TEMP1*(MESH(J,4)-MESH(J-1,4)) 
MESH(I+1,2)  =  MESH(J-1.5)  +  TEMP1*(MESH( J,5) -MESH( J-1 ,5 ) ) 
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ENDIF 

60  CONTINUE 

ENDIF 
C 
C 

C  BISECTION  NODE  INITIALIZATION  (X,Y  COORD) 

C 

WRITE(6,*)  'BISECTION  NODE  #  =  ’.BIND 

IF( (METHOD.  EQ.  1).  OR.  (METHOD.  EQ.  3))  THEN 
DO  70,  1=1,  BIND 

TEMPI  =  I/FLOAT(BIND  +1) 

MESH(I+PERND,1)  =  MESH(1,1)  +  TEMPl*(MESH(BIND+2, 1)  - 

CMESH(1,1)) 

MESH(I+PERND,2)  =  MESH(1,2)  +  TEMPl*(MESH(BIND+2 , 2)  - 

CMESH(1,2)) 

70  CONTINUE 

C 

ELSE 

DO  80,  1=2,  PERND+1 

MESH(PERND+I-1,1)  =  MESH( PERND+2 - 1 , 1 )  +  0.  5*(MESH( 1 , 1) 
C  MESH(PERND+2-I,l)) 

MESH(PERND+I-1,2)  =  MESH(PERND+2-I ,2)  +  j. 5*(MESH( I , 2) 
C  MESH(PERND+2-I,2)) 

80  CONTINUE 

C 

C 

ENDIF 

C 

IF((METHOD.EQ.  5).  OR.  (METHOD.  EQ.  6))  THEN 
DO  90,  M  =  1,  BIND 

MESH(M,4)  =  MESH(PERND+M,1) 

MESH(M,5)  =  MESH(PERND+M,2) 

90  CONTINUE 

DISTB  =  0.0 
MESH(0,3)  =  0.  0 
DO  100,  L  =  1,  BIND+1 
IF(L.EQ.  1)  THEN 

DIST  =  SQRT((MESH(1,1)  -  MESH( 1,4))**2  + 
C(MESH(1,2)  -  MESH(1,5))**2) 

ELSEIF(L.  EQ. BIND+l)  THEN 

DIST  =  SQRT((MESH(BIND,4)  -  MESH( BIND+2 , 1) ) 

C**2  +  (MESH(BIND,5)  -  MESH(BIND+2,2) )**2) 

ELSE 

DIST  =  SQRT((MESH(L-1,4)  -  MESH(L,4)) 

C**2  +  (MESH(L-1,5)  -  MESH(L,5))**2) 

ENDIF 

DISTB  =  DISTB  +  DIST 
MESH(L,3)  =  DISTB 
100  CONTINUE 

DZ  =  DISTB/(BIND  +  1.0) 

WRITE(6,*)  'DZ  SPACING  =  ' ,  DZ 

C 

DO  120,  I  =  1,  BIND 
TEMP  =  DZ*I 
J  =  1 

110  IF(TEMP.  GT. MESH(J,3))  THEN 


86 


o  o  o  n  n  n 


J  =  J  +  1 
GOTO  110 

ENDIF 

IFCJ.EQ.  1)  THEN 

TEMPI  =  TEMP/ (SQRT((MESK( 1,1)  -  MESH( 1 ,4) )**2  + 
C  (MESH(1,2)  -  MESH(1,5))**2)) 

MESH(PERND+I,1)  =  MESH(1,1)  +  TEMP 1*(MESK( 1,4) 

C  -  MESH(1,1)) 

MESH(PERND+I,2)  =  MESH(1,2)  +  TEMP1*(MESH( 1 ,5) 

C  -  MESH(1,2)) 

ELSE 

TEMPI  =  (TEMP  -  MESH(J-1,3))/(SQRT((MESH(J,4)  > 

C  MESH(J-1,4))**2  +  (MESH(J,5)  -  MESH( J-l,5))**2)) 

MESH(PERND+I,1)  =  MESH(J-1,4)  +  TEMP1*(MESH( J,4) 
C  -  MESH(J-1,4)) 

MESH(PERND+I,2)  =  MESH(J-1,5)  +  TEMP1*(MESH( J,5) 
C  -  MESH(J-1,5)) 

ENDIF 

120  CONTINUE 

ENDIF 
RETURN 
END 


SUBROUTINE  NORMAL 

THIS  ROUTINE  COMPUTES  THE  X  AND  Y  COMPONENTS 
OF  THE  OUTWARD  UNIT  NORMAL  AT  EACH  SURFACE  POINT. 

REAL  DR,  DZ,  DL,  MESH(0; 200,5) 

INTEGER  I.  PERND,  BIND 
COMMON/BLKl/MESH 
C0MM0N/BLK2/PERND,BIND 
DO  10,  1=1,  PERND 
IF(I.EQ.  1)  THEN 

DR  =  MESH(2,1)  -  MESH( PERND, 1) 

DZ  =  MESH(2,2)  -  MESH(PERND,2) 

DL  =  SQRT(DR*DR+DZ*DZ) 

MESH(1,3)=-DZ/DL 
MESH(1,4)=DR/DL 
ELSEIFC I.  EQ.  PERND)  THEN 

DR  =  MESH(1,1)  -  MESH( PERND -1,1) 

DZ  =  MESH(1,2)  -  MESH( PERND -1,2) 

DL  =  SQRT(DR*DR+DZ*DZ) 

MESH( PERND , 3 )=-DZ/DL 
MESH(PERND,4)=DR/DL 

ELSE 

DR  =  MESH(I+1,1)  -  MESH(I-1,1) 

DZ  =  MESH(I+1,2)-  MESH(I-1,2) 

DL  =  SQRT(DR*DR+DZ*DZ) 

MESH(I,3)=-DZ/DL 
MESH(I,4)=DR/DL 

ENDIF 
10  CONTINUE 
RETURN 
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IHD 


SUBROUT I.VE  NOOStT(  «ET«00 . SABC . SIMX . UMt ) 

THIS  ROUTINE  USES  A  TVi’O-SVEEP  TECHNIQUE  Tv»  CCWPITE  THE 
NUMBER  OF  NODES  ALONG  EACH  N)^C  ROH,  I  A  fOLSH 
ORIENTATION  ATTRIBUTE  IS  ALSO  SET 

SET  Z-AXIS  SPACING  AND  LSDPOINT  NODES 

REAL  DZ.  ZZ.  HESKIO:  200.5).  0.  DIST.  OISTB 
INTEGER  I.  J.  NO0ES(20O).  MORI 200).  PERNS.  SOLD.  NNCV.  BIND 
INTEGER  L.  METHOD.  NABC( 100.3).  NBMX.  UNR 
CHARACTER*!  CHAR2 .  CHAR3.  CHAR4 
COMMON/ BUCl /MESH 
COMMON/ B  LK2 / PERND .BIND 
COMMON/ BIJC3/ NODES 
COMMON/ BLK6/M0R 
COMMON / B  LK9 / CHAR2 . CHAR  3 . CKAR4 

UNK  ■  0 
NBMX  ■  -99 

IF(( METHOD.  EQ.  1).  OR.  (METHOD. EQ  3))  THEN 

DZ  ■  (SQRT((MESH(  1.2)  •  MESH(  BlSt>*2 . 2 )  )►•:  ♦ 

C  (MESHd.l)  •  MESK(BINI>*2.)))**2))/(81NO  •  1  0) 

ELSE 

DISTB  -0.0 
DO  10.  L  -  1.  BINIHl 
IF(L.  EQ.  1)  THEN 

DIST  ■  SQRT(  ( MESH(  1 . 1 )  •  MESH(  PERM)*! ,  1)  )**2  * 
C(MESH(l,2)  -  MESH(PER.Nt>*l,2))**2) 

ELSEIF(L.  EQ.  BIND*!)  THEN 

DIST  -  SQRT((MESH(PER.VD*BINO.!)  •  MESH(  £  JND*:  , !  ;  y 
C**Z  +  (MESH(PERND*BIN-0,2)  -  MESH(BIND*2,2)  )**2) 

ELSE 

DIST  «  SQRT({MESH(PER.S*D*I-1.1)  •  MESH(  PERM)*  1  . !  ) ) 
C**2  +  (MESH(PERND+I-1,2)  -  MESH( PERM)*1 , 2 ) )**2 ) 

END  IF 

DISTB  «  DISTB  ♦  DIST 
10  CONTINUE 

DZ  =  DISTB/(BIND  +  1.0) 

ENDIF 

WRITE(6,*)  'BISECTION  SPACING  «  ’ ,  DZ 
NODES(l)  =  2 
N0DES(2)  =  3 
N0DES(BIND+2)  =  2 
N0DES(PERND+1)  =  2 
NODES  (PERM))  =  3 
C  PERFORMING  FORWARD-SWEEP 
DO  20  1=3,  BIND+1 

NOLD  =  NODES(I-l) 

D  =  SQRT(( MESHd.l)  -  MESH( I+PERND-1 , 1) )**2  + 

C  (MESH(I,2)-MESHd+PERND-l,2))**2) 

NNEW  =  INT(0.  1  +  D/DZ)  +  2 
NODES(I)  =  NNEW 
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IF  (NNEW.  GT.  NOLIHl)  NODES(I)  =  NOLD  +  1 
IF  (KNEW.  LT.  NOLD- 1)  NODES(I)  =  NOLD  -  1 
IF  (NODES(I).  LT.  3)  NODES(I)  =  3 
20  CONTINUE 
C 

J  •  PERND  +  2 

DO  30.  I  ■  PERND-1,  BIND+3,  -1 
NOLD  ■  NODES(I+l) 

D  ■  SQRT((NESH(I,1)  -  MESH( J,l))**2  + 

C  (MESH(I.2)-MESH(J,2))**2) 

NNEV  «  INT(0.  1  +  D/D2)  +  2 
NODES(I)  ■  NNEW 

IF  (NNEW.  GT.  NOLD+1)  NODES(I)  =  NOLD  +  1 
IF  (NNEW. LT. NOLD- 1)  NODES(I)  =  NOLD  -  1 
IF  (NODES(I).  LT.  3)  NODES(I)  =  3 
J  -  J  +  1 
30  CONTINUE 
C 

C  BACK-SWEEP  TO  RESET  LAST  NODES  IF  NEEDED 

C 

I  -  BIND  ♦  2 

40  I  -  I  -  1 

IF  (I.  EQ.  2)  GO  TO  50 

IF  (NODES(I).  LE.NODES(lTi)-rl)  GO  iO  60 
NODES(I)  ■  N0DES(I+1)  +  1 
GO  TO  40 
50  CONTINUE 

WRITE(6,*)  '  PROGRAM  ABORTED  IN  NODSET  RIGHT  SIDE  BACKSWEEP  ’ 
STOP 

60  CONTINUE 

C 

I  ■  BIND  +  2 

70  I  »  I  +  1 

IF  (I.EQ.  PERND)  GO  TO  80 
IF  (NODES(I).  LE.N0DES(I-1)+1)  GO  TO  90 
NODES! I)  -  NODES! I-l)  +  1 
GO  TO  70 
80  CONTINUE 

WRITE! 6,*)  ’  PROGRAM  ABORTED  IN  NODSET  LEFT  HALF  BACKSWEEP  ' 
STOP 

90  CONTINUE 

C 

C  FORWARD  SWEEP  TO  LOAD  MESH  ORIENTATION  ARRAY,  MOR 

C 

MOR!l)  »  0 
MOR!BIND+2)  -  0 
M0R!PERNI>+1)  »  0 
C 

DO  100,  I»2,  BIND+1 

IF!NODES!I+l).GT. NODES!!))  THEN 
MOR! I)  *  0 

ELSEIF(NODES(I+l).  LT.  NODES(I))  THEN 
MOR! I)  =  1 

ELSE 

IF!NODES!I+2).GT.NODES!I))  THEN 
MOR! I)  =  0 
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ELSEIF(NODES( 1+2) .  LT.  NODES( I ) )  THEN 
MOR(I)  =  1 

FT  ^F 

MOR(I)  =  MOR(I-l) 

ENDIF 

ENDIF 

100  CONTINUE 
C 

DO  110,  I  =  PERND,  BIND+3,  -1 

IF(NODES(I-l).GT.NODES(I))  THEN 
MOR(I)  =  0 

ELSE IF( NODES ( I - 1 ) .  LT.  NODES ( I ) )  THEN 
MOR(I)  =  1 

ELSE 

IF(NODES(I-2).GT.NODES(I))  THEN 
MOR(I)  =  0 

ELSEIF(NODES(I-2).  LT. NODES( I))  THEN 
MOR(I)  =  1 

ELSE 

MOR(I)  =  MOR(I+l) 

ENDIF 

ENDIF 

110  CONTINUE 
C  LOAD  NABC  ARRAY 

DO  120,  I  =  1,  BIND+2 
IF(I.EQ. 1)  THEN 

NABC(1,1)  =  0 
NABC(1,2)  =  1 
NABC(1,3)  =  3 
ELSEIFCI.  LE.  BIND)  THEN 

NABC(I,1)  =  NABC( 1-1,2) 

NABC(I,2)  =  NABC(I-1,3) 

NABC(I,3)  =  NODESC PERND- I+l)  +  NODES(I+l)  -  3 
ELSEIFCI. EQ.  BIND+1)  THEN 

NABC(I,1)  =  NABC(I-1,2) 

NABC(I,2)  =  NABC(I-1,3) 

NABC(I,3)  =  1 

ELSE 

NABC(I,1)  =  NABC(I-1,2) 

NABC(I,2)  =  NABCC 1-1,3) 

NABC(I,3)  =  0 

ENDIF 

120  CONTINUE 
C 

DO  130,  1=1,  BIND+2 

IF((CHAR2.EQ.  'l').0R.  (CHAR2.  EQ.  '  i' )  )  THEN 

WRITE(12,*)  I, NABCC 1,1), NABCC I, 2), NABCC I, 3) 

ENDIF 

UNK  =  UNK  +  NABCC  1, 2) 

IFCNABCCI,2).GE.NBMX)  THEN 
NBMX  =  NABCC I, 2) 

ENDIF 

130  CONTINUE 
C 

WRITEC*,*)  'MAXIMUM  UNKNOWN  WIDTH  =  ’ ,NBMX 
WRITEC*,*)  'TOTAL  #  OF  UNKNOWNS  =  ' ,UNK 
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WITSC*.*)  •  ’ 

RETUtN 

IXD 


SUBROVriNE  SOrmid.LEL.LAND) 

THIS  SUlRCKTriNI  CEKLRATTS  A  HESH  ROW  FOR  THE  INPUT  ROW  I. 

LOADIW  or  THE  LOCAL  NODE -ELEMENT  CWWECTION  MATRICES  LND  AND 
SDL  FOR  ELEMENTS  SETVEEN  I  AND  !♦!  NODE  ROWS.  REFERENCE  IS  TO 
THE  LEFT  SIDE  OF  THE  1th  ROW  OR  VECTOR. 

INTEGER  N0DES(200).  MOR(200).  LKD( 0: 200.3) ,  NDL(200,4),  NCT(200) 

INTEGER  I.  PER.ND.  BIND,  LEL,  LAND,  J.  K.  LL,  JJ,  NN,  KK,  N,  Nl,  N2 

INTEGER  NXyMX.  NOSIL.  NOS2L.  KDSIR.  NDS2R.  LMX 

COMMON/ BUC2  /  PERND .  B 1 ND 

COMMON/ BUC3/SODES 

COMMON/ BLIC4/ LND .  SDL .  NCT 

COMMON*/ BLKB/MOR 

IF(  I.  £Q.  BINIH2)  THEN 

WRITE(6.'*)  ‘ERRORED  OUT  IN  SUBROUTINE  SORTER,  YOU  ATTEMPTED' 
WRITE! 6.*)  '  TO  CALL  SORTER  WITH  I  -  BIND  +  2!  ' 

WRITE! A.*)  '  THIS  ROW  HAS  NO  ELEMENTS' 

RETURN 

ENDIF 

DO  20,  J  -  0,  200 

DO  10.  K  ■  1.  3 
LND!J.IC)  ■  0 
0  CONTINUE 

0  C0^f^ISUE 
NDMX  ■  200 

N-DSIL  -  NODES!PERNrH2-I) 

NDS2L  -  NODES! PERNIHl-I) 

SDSIR  «  NODES! I) 

NDS2R  -  NODESII+l) 

LMX  ■  NDSIL  ♦  NDSIR  ♦  NDS2L  ♦  NDS2R  -  2 


LEFTSIDE  OF  THE  BISECTION  SEGMENT 
TOP  ROW 

IF! I. EQ. 1)  THEN 
LND! 1,1)  »  4 
LND!1,2)  ■  1 
LND!1,3)  «  3 
LND! 2,1)  «  1 
LND!2,2)  «  4 
LND!2.3)  «  2 
LND(3,1)  »  5 
LND!3,2)  «  2 
LND!3,3)  a  4 
LND!4,1)  a  5 
LND!4,2)  a  2 
LND!4,3)  a  6 
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LND(5,1) 
LND(5,2) 
LND(5,3) 
LND(6,1) 
LND(6,2) 
LND(6,3) 
LEL  =  6 
LAND  =  7 


BOTTOM  ROW 

ELSEIF(I.EQ.  BIND+1)  THEN 


LND(1,1) 

LND(1,2) 

LND(1,3) 

LND(2,1) 

LND(2,2) 

LND(2,3) 

LND(3,1) 

LND(3,2) 

LND(3,3) 

LND(4,1) 

LND(4,2) 

LND(4,3) 

LND(5,1) 

LND(5,2) 

LND(5,3) 

r.ND(6,l) 

LND(6.2) 

LND(6,3) 

LEL  =  6 

LAND  =  7 


2 

6 

1 

6 

2 

7 

3 

7 

2 

3 
7 

4 
6 
4 
7 

4 
6 

5 


30 

40 

C 


EQUAL  NODE  NUMBERS 

ELSEIFCNDSIL.  EQ.NDS2L)  THEN 
FOR  MOR  =  0  (LH  ORIENTATION) 

IF(M0R(PERND+2-I).EQ.  0)  THEN 
LND(1,1)  =  1 

LND(I,2)  =  NDSIL  +  NDSIR 
LND(1,3)  =  2 

LND(2,1)  =  NDSIL  +  NDSIR  +  1 
LND(2,2)  =  2 

LND(2,3)  =  NDSIL  +  NDSIR 
DO  40,  N  =  1,  NDSIL-2 
N1  =  2*N  +  1 
N2  =  N1  +  1 
DO  30,  K  =  1,  3 

LND(N1,K)  =  LND(1,K)  +  N 
LND(N2,K)  =  LND(2,K)  +  N 
CONTINUE 
CONTINUE 

1  (RH  ORIENTATION) 


FOR  MOR  = 
ELSE 


LND(1,1)  =  NDSIL  +  NDSIR 
LND(1,2)  =  1 
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LKD(1.3)  •  KDSIL  ♦  SDSIR  ♦  I 
LND(2.l)  •  2 

LVD(2.2)  •  KOSIL  ♦  KDSIR  ♦  I 
LKO(2.3)  •  I 
DO  60.  N  »  I,  W)SIL*2 
Nl  ■  2*N  ♦  I 
N2  ■  Nl  ♦  I 
DO  50.  K  •  I,  3 

LND(NI.K)  ■  LKD(l.K)  ♦  N 
LND(N2.K)  ■  LKD(2.K)  ♦  N 
50  COffTINLE 

60  CO^frI^^X 

ENDIF 
C 

C  LEmUND  MFSH  ORIENTATION 
C 

ELSElF(«OR(PERND+2-I).  EQ.  0)  THEN 
LND(l.l)  -  NOSIL  ♦  NDSIR  1 
LND(1.2)  -  I 

LND( 1.3)  -  NDSIL  ♦  NDSIR 
LND(0,l)  -  0 

LVD(0,2)  ■  SDSIL  ♦  NDSIR 
LND(0.3)  -  I 
DO  80.  N  ■  1.  NDSIL- 1 
Nl  ■  2*N 
N2  -  Nl  ♦  1 
DO  70.  K  -  1.  3 

LND(Nl.K)  ■  LND(0,K)  ♦  N 
LN0(N2.K)  -  LSDd.K)  +  N 
70  COST  I  NIX 

80  CONTIN-UE 
C 

C  RIGHTHAND  MESH  ORIENTATION 

C 

ELSE 

LNDd.l)  »  2 

LNDd,2)  »  NDSIL  ♦  NDSIR 
LNDd,3)  «  1 

LND(0,1)  «  NDSIL  +  NDSIR  •  1 
LND(0,2)  -  I 

LND(0,3)  =  NDSIL  +  NDSIR 
DO  100,  N  »  1,  NDSlL-2 
Nl  *  2*N 
N2  *  Nl  +  1 
DO  90,  K  ■  1,  3 

LND(N1,K)  «  LND(0,K)  +  N 
LND(N2,K)  »  LNDCl.K)  +  N 
90  CONTINDX 

100  CONTINUE 
ENDIF 
C 

IF((I.EQ.  D.OR.  (I.EQ.  BIND+D)  THEN 
GOTO  230 

ENDIF 

C 
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LEL  =  N2 


C 

C 

C  RIGHTSIDE  OF  THE  BISECTION  SEGMENT 

C 

C 

C  EQUAL  NODE  NUMBERS 

C 

IF(NDS1R.  EQ.NDS2R)  THEN 
C  MOR  =  0  (LH  ORIENTATION) 

IF(MOR(I).EQ.  0)  THEN 

LND(LEL+1,1)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 
LND(LEL+1,2)  =  NDSIL 
LND(LEL+1,3)  =  NDSIL  +  NDSIR  +  NDS2L 
LND(LEL+2,1)  =  NDSIL  +  1 
LND(LEL+2,2)  =  NDSIL  +  NDSIR  +  NDS2L 
LND(LEL+2,3)  =  NDSIL 
DO  120,  N  =  1,  NDSlR-2 
N1  =  2*N  +  1  +  LEL 
N2  =  N1  +  1 

DO  110,  K  =  1,  3 

LND(N1,K)  =  LND(LEL+1,K)  +  N 
LND(N2,K)  =  LND(LEL+2,K)  +  N 
110  CONTI.NUE 

120  CONTINUE 

LEL  =  N2 
LAND  =  LMX 

C  MOR  =  1  (RH  ORIENTATION) 

ELSE 

LND(LEL+1,1)  =  NDSIL 

LND(LEL+1,2)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 
LND(LEL+1,3)  =  NDSIL  +  1 
LND(LEL+2,1)  =  NDSIL  +  NDSIR  +  NDS2L 
LND(LEL+2,2)  =  NDSIL  +  1 
LND(LEL+2,3)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 
DO  140,  N  =  1,  NDSlR-2 
N1  =  2*N  +  1  +  LEL 
N2  =  N1  +  1 

DO  130,  K  =  1,  3 

LND(N1,K)  =  LND(LEL+1,K)  +  N 
LND(N2,K)  =  LND(LEL+2,K)  +  N 
130  CONTINUE 

140  CONTINUE 

ENDIF 
LEL  =  N2 
LAND  =  LMX 
C 

C  LEFT  HAND  MESH  ORIENTATION 
C 

ELSEIF(MOR(I).EQ.  0)  THEN 

LND(LEL+1,1)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 

LND(LEL+1,2)  =  NDSIL 

LND(LEL+1,3)  =  NDSIL  +  NDSIR  +  NDS2L 

LND(LEL+2,1)  =  NDSIL  +  1 

LND(LEL+2,2)  =  NDSIL  +  NDSIR  +  NDS2L 
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LND(LEL+2,3)  =  NDSIL 

DO  160,  N  =  1,  NDSlR-1 
N1  =  2*N  +  LEL  +  1 
DO  150,  K  =  1,  3 

LND(N1,K)  =  LND(LEL+1,K)  +  N 
CONTINUE 
CONTINUE 

DO  180,  N  =  1,  NDSlR-2 
N2  =  2*N  +  LEL  +  2 
DO  170,  K  =  1,  3 

LND(N2,K)  =  LND(LEL+2,K)  +  N 
CONTINUE 
CONTINUE 

LEL  =  N1 
LAND  =  LMX 

RIGHT  HAND  MESH  ORIENTATION 


ELSE 


LND(LEL+1,1)  =  NDSIL 

LND(LEL+1,2)  =  NDSIL  +  NDSIR  +  NDS2L  - 

LND(LEL+1,3)  =  NDSIL  +  1 

LND(LEL+2,1)  =  NDSIL  +  NDSIR  +  NDS2L 

LND(LEL+2,2)  =  NDSIL  +  1 

LND(LEL+2,3)  =  NDSIL  +  NDSIR  +  NDS2L  - 

DO  200,  N  =  1,  NDSlR-2 
N1  =  2*N  +  LEL  +  1 
DO  190,  K  =  1,  3 

LND(N1,K)  =  LND(LEL+1,K)  +  N 
CONTINUE 
CONTINUE 

DO  220,  N  =  1,  NDSlR-3 
N2  =  2*N  +  LEL  +  2 
DO  210,  K  =  1,  3 

LND(N2,K)  =  LND(LEL+2,K)  +  N 
CONTINUE 
CONTINUE 

LEL  =  N1 
LAND  =  LMX 


ENDIF 


ZEROING  NDL  AND  NCT  PRIOR  TO  FILL 

DO  250,  N  =  1,  NDMX 
NCT(N)  =  0 
DO  240,  K  =  1,  4 


NDUN.K)  »  0 
240  COVTINLT 

250  CONTINUE 
C 

C  SCANNING  LND  ARJIAV  TO  LOAD  M)L 
C 

DO  270  LL  -  1.  LEL 

DO  260  JJ  ■  1.  3 

NN  •  LND(LL.JJ) 

NCT(NN)  •  NCT(NN)  ♦  I 
KK  ■  NCT(NN) 

NDUNN.KX)  •  LL 
260  CONTINUE 

270  CONTINUE 
C 

END 

C 

c 

SUBROUTINE  FINDER( I . LAND .DIST) 

C 

C  THIS  SUBROUTINE  DETERMINES  THE  (X.Y)  COORDINATES  OF  EACH 
C  NODE  IN  THE  CALLING  Ith  ROV  OR  VECTOR  MESH. 

C 

INTEGER  I.  J,  LAND.  PERND,  BIND.  NODES! 200) 

REAL  MESH(0:  200.5).  NOP( 200.2).  DIST 

COMMON/BLKl/MESH 

COMMON/  B  LK2  /  PERNT) .  B I ND 

COMMON/ BLK3/NODES 

C0MM0N/BLK5/NDP 

iFd.EQ.  I)  THEN 

NDP(l.l)  ■  MESH(l.l)+MESH(1.3)*DIST 
NDP(1,2)  -  MESH(  1.2)+MESH(  1.4)*DIST 
NDP(2.l)  =  MESH! 1,1) 

N’DP(2,2)  -  MESHC  1,2) 

NDP(3,1)  «  MESH( PERND. 1)+MESH( PERND, 3 )*DI ST 
NDP(3.2)  -  MESH(PER.ND,2)-*’MESH(PERNT),4)*D1ST 
NDP(4,1)  -  MESH(PERND,1) 

NDP(4,2)  «  MESH(PERNT).2) 

NDP(5,1)  «  MESH(PERNT>+1,1) 

NDP(5,2)  «  MESH(PERND+1.2) 

NDP(6,1)  »  MESH(2,1) 

NDP(6,2)  «  MESH(2,2) 

NDP(7,1)  «  MESH(2,1)+MESH(2,3)*DIST 
NDP(7,2)  *  MESH{2,2)+MESH(2,4)*D1ST 
ELSEIFCI.EQ.  BIND+1)  THEN 

NDP(1,1)  »  MESH(BINI>+3,1)+MESH(BIND+3,3)*DIST 
NDP(1,2)  =  MESH(BIND+3,2)+MESH(BIND+3,4)*DIST 
NDP(2,1)  «  MESHC BINIH3,1) 

NDP(2,2)  »  MESH(BIND+3,2) 

NDP(3,1)  »  MESHC PERND+BIND.l) 

NDP(3,2)  «  MESH(PERND+BIN'D,2) 

NDP(4,1)  =  MESHC BIND+l.l) 

NDPC4,2)  =  MESHC BIND+1, 2) 

NDPC5,1)  =  MESHC BIND+1, 1)+MESHC BIND+1, 3)*DIST 
NDPC5,2)  =  MESHCBIND+1,2)+MESHCBIND+1,4)*DIST 
NDPC6,1)  =  MESHC BIND+2,1)+MESHCBIND+2,3)*DIST 
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NDP(6,2)  =  MESH(BIND+2,2)+MESH(BIND+2,4)*DIST 
NDP(7,1)  =  MESH(BIND+2.1) 

NDP(7,2)  =  MESH(BIND+2.2) 

ELSEIF( I.  EQ. BIND+2)  THEN 

WRITE(6,*)  '  ERRORED  OUT  IN  SUBROUTINE  FINDER,  YOU  ATTEMPTED* 
WRITE(6,*)  '  TO  CALL  FINDER  WITH  I  =  BIND  +  2!  ' 

WRITE(6,*)  '  THIS  ROW  HAS  NO  ELEMENTS  AND  NO  COORDINATES' 

ELSE 

C  NODE  1 

NDP(1,1)  =  MESH(PERND-I+2,1)+MESH(PERND-I+2,3)*DIST 
NDP(1,2)  =  MESH(PERND-I+2,2)+MESH(PERND-I+2,4)*DlST 
C  NODE  2  TO  THE  BISECTION  SEGMENT 

DO  10,  J  =  2,  N0DES(PERND“I+2) 

NDP( J , 1 ) =MESH( PERND - 1+2 , 1 ) - ( J - 2 ) *( MESH( PERND - 1+2 , 1 ) - 
C  MESH(PERND+I-l,l))/(N0DES(PERND-I+2)-2) 

NDP( J , 2 )=MESH( PERND- 1+2 , 2 ) - ( J- 2 )*( MESH( PERND- 1+2 , 2 ) - 
C  MESH( PERND+I - 1 , 2 ) ) / ( NODES( PERND- 1+2 ) - 2 ) 

10  CONTINUE 

C  BISECTION  SEGMENT  TO  THE  RIGHTSIDE  SURFACE 

DO  20,  J  =  3,  NODES(I) 

NDP( NODES( PERND-I+2 )+J-2 , 1 )=MESH( PERND+I - 1 , 1 )+( J-2 )* 

C  (MESH( I , 1) -MESH(PERND+I-1, 1))/(N0DES( I)-2) 

NDP(N0DES(PERND-I+2)+J-2,2)=MESH(PERND+I-l,2)+(J-2)* 

C  (MESH(I,2)-MESH(PERND+I-l,2))/(N0DES(I)-2) 

20  CONTINUE 

C  Ith  ROWS  LAST  NODE 

NDP(NODES(PERND-I+2)+NODES(I)-l,l)=MESH(I,l)+MESH(I,3)*DIST 
NDP( NODES ( PERND-I+2 ) +NODES( I ) - 1 , 2 ) =MESH ( I , 2 ) +MESH ( I , 4 ) *DI ST 
C  I+lth  ROWS  FIRST  NODE 

NDP( NODESC  PERND- 1+2 )+NODES( I ) , 1)*MESH( PERND- I+l , 1 )+MESH 
C  (PERND-I+1,3)*DIST 

NDPC  NODESC  PERND- 1+2 )+NODES( I ) , 2 )=MESH( PERND- I+l , 2 )+MESH 
C  (PERND-I+1,4)*DIST 

C  I+ITH  ROW  (NODE  2)  TO  THE  BISECTION  SEGMENT 

DO  30,  J  =  2,  NODESC PERND-I+1) 

NDPC  J+NODES( PERND-I+2)+NODES( I ) - 1 , 1)=MESH( PERND-i+1 , 

C  1 ) - ( J-2 )*( MESHC  PERND-I+1 , 1 ) -MESH( PERND+I , 1) ) / 

C  (NODESC PERND- I+l) -2) 

NDPC  J+NODE  S ( PERND - 1+2 ) +NODES ( I ) - 1 , 2 ) =MESH( PERND - 1+ 1 , 

C  2)-(J-2)*(MESH(PERND-I+l,2)-MESH(PERND+I,2))/ 

C  (NODESC PERND-I+1) -2) 

30  CONTINUE 

C  I+lth  ROW  BISECTION  SEGMENT  TO  THE  RIGHTSIDE  SURFACE 

DO  40,  J  =  3,  NODESC I+l) 

NDPC LAND-NODES(I+l)+J-l,l)=MESH(PERND+I,l)+( J-2)* 

C  ( MESHC 1+1,1) -MESHC  PERND+I , 1 ) )/( NODESC I+l ) -2 ) 

NDP(LAND-NODES(I+l)+J-l,2)=MESH(' PERND+I,  2)+(  J-2)* 

C  ( MESHC I+l , 2 ) -MESHC  PERND+I , 2 ) ) / ( NODESC 1+ 1 ) -2 ) 

40  CONTINUE 

C  LAST  NODE 

NDP(LAND,1)  =  MESHC I+l, 1)+MESH( I+l, 3)*DIST 
NDP(LAND,2)  =  MESHC I+l ,2)+MESH( I+l ,4)*DIST 
C 

ENDIF 

RETURN 
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no  o  t-  o  n  o  n  o  o  o  o  n  o  o  o  o  n  o  no 


END 


SUBROUTINE  VARINT( J , F , ALPHA , BETA , AREA , LND) 

GENERATING  VARIATIONAL  FINITE  ELEMENT  AREA  INTEGRATIONS  OF  THE 
LINEAR  BASIS  FUNCTION  LAGRANGIAN  FOR  THE  HELMHOLTZ  EQUATION. 
THESE  ARE  RETURNED  IN  F(3,3).  X(3)  AND  Y(3)  ARE  THE  WAVENUMBER 
NORMALIZED  CARTESIAN  COORDINATES  OF  THE  TRIANGLE  VERTICES. 

X  =  Ko*x,  Y  =  Ko*y  ALPHA  AND  BETA  ARE  COMPLEX 
MATERIAL  PARAMETERS  WITHIN  THE  ELEMENT. 

FOR  TM  INCIDENCE:  ALPHA  =  1/ur;  BETA  =  er 
FOR  TE  INCIDENCE:  ALPHA  =  1/er;  BETA  =  ur 

OUTPUTS  :  F(3,3)  -  FINITE  ELEMENT  AREA  INTEGRATION 

AREA  -  AREA  OFF  A  TRIANGLE 

COMPLEX  ALPHA,  BETA,  B12,  F(3,3) 

REAL  NDP(200,2),  X(3),  Y(3),  T(3,3),  AREA,  DET 
INTEGER  L,  K,  J,  LND(0: 200,3) 

COMMON/BLK5/NDP 

X(l)  =  NDP(LND(J,1),1) 

X(2)  =  NDP(LND(J,2),1) 

X(3)  =  NDP(LND(J,3),1) 

Y(l)  =  NDP(LND(J,1),2) 

Y(2)  =  NDP(LND(J,2),2) 

Y(3)  =  NDP(LND(J,3),2) 

DET  =  X(2)*Y(3)  +  X(3)*Y(1)  +  X(1)*Y(2)  -  X(3)*Y(2)  - 
CX(1)*Y(3)  -  X(2)*Y(1) 

AREA  =  ABS(0.5*DET) 

B12  =  BETA/12. 

T(l,l)  =  (Y(2)  -  Y(3))/DET 

T(l,2)  =  (Y(3)  -  Y(1))/DET 

T(l,3)  =  (Y(l)  -  Y(2))/DET 

T(2,l)  =  (X(3)  -  X(2))/DET 

T(2,2)  =  (X(l)  -  X(3))/DET 

T(2,3)  =  (X(2)  -  X(1))/DET 

T(3,l)  =  (X(2)*Y(3)  -  X(3)*Y(2))/DET 

T(3,2)  =  (X(3)*Y(1)  -  X(1)*Y(3))/DET 

T(3,3)  =  (X(1)*Y(2)  -  X(2)*Y(1))/DET 

DO  10,  K  =  1,  3 

DO  10,  L  =  1,  3 

F(K,L)  =  ALPHA*(T(1,K)*T(1,L)  +  T( 2 ,K)*T( 2 ,L) )  -  B12 
IF(K.  EQ.  L)  F(K,L)  =  F(K,L)  -  B12 
0  F(K,L)  =  AREA*F(K,L) 

RETURN 

END 


SUBROUTINE  LOADER( BCOND , OFFSET, ALPHA , BETA , NABC , IMX ,NBMX , EO , SURBC , 
CCHARl , LINE , MODE , XORIGIN , YORIGIN , SVEC , CHARS ) 

C 

COMPLEX  A(50,50),B(50,50),C(50,50),P(50,100) 
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COMPLEX  F(3,3),FROV(  100 . 3 . 3)  .BCOS'D(  100)  .LISEl  50) 

COMPLEX  ALPHA,  BETA,  DETERM,  SURBC(IOO),  SVECOO.lOO) 

REAL  OFFSET, MINAREA,MAXAREA, AREA. RATIO, EO.KJ.XCRIGIS.YOSIGIN 
REAL  UBCONO 

INTEGER  I , J , JD . K , L . NDTOP . NDBOT . NDTOT , NOD . KND . KD(  3 ) . LEL . LAND 
INTEGER  LND(0;  :’''0.3).  NDL(200.4).  NCTCOO).  PERKD,  BIND.  JJ 
INTEGER  NODES! 200),  HISROW.  MINEL,  MAXROV,  HAXEL,  TCALL,  NL 
INTEGER  N , M . NBMX . NABC( 100.3), INORM . I MX . MODE 
CHARACTER*!  CHARI,  CHAR2 .  CHAR3,  CHAK4.  CHARS 
C 

COMMON/ BLK2 / PERND .BIND 
C0MM0N/BLK3/N0DES 
COMMON  /  B  LK4  /  LND ,  NT)  L .  SCT 

COMMON  /  B  LK  7  /  M I N  ARE  A .  M I NROW .  M I NE  L .  MAXAREA ,  MAXROb- .  MAXE  L .  AREA 
COMMON/BLK8/A.B.C.P 
COMMON /BLK9/ CHAR 2.  CHAR3.  CHAR4 
C 

UBCOND  -  1. 0 
INORM  •  0 
IMX  -  BINT)  ♦  2 
I  »  1 

TCALL  ■  BIND  +  i 
C 

IF((CHAR3.  EQ.  ’D’  ).  OR.  (CHAR3,  EQ.  ‘d’  ))  THEN 
WRITE! 20. 1030) 

WRITE! 20. 1040) 

WRITE!20. 1050) 

WHITE! 20. 1060) 

WRITE!20,1070) 

END  IF 

WHITE!*, 1000)  I, TCALL 
C 

CALL  SORTER!!. LEL. LAND) 

CALL  FINDER! I .LAND, OFFSET) 

IF!. NOT.  !! CHARS. EQ.  'M' ).  OR. ! CHARS. EQ. ’»’)))  THEN 
CALL  BNDC! I .BCOND.EO.SIHBC, CHARI .ALPHA, BETA, LISE, MODE, XORIGIS, 
CY0RIGIN,CHAR4) 

END  IF 

CALL  FILL! I.LEL,FROW, ALPHA. BETA) 

IF!.NOT.  !!CHAR5.EQ.  ’M’ ).  OR.  !CHAR5.  EQ.  ’m’)))  THEN- 
CALL  ZERO 
END  IF 

C  ESTABLISH  THE  NUMBER  OF  NODES  IN  THE  I  ■  1  AND  I  *  2  ROWS 
NDTOP  »  2 
NDBOT  =  5 
NDTOT  =  7 
C 

C  B,  C  &  P  FOR  I  »  1  !TOP) 

C 

J  =  1 

DO  70,  JD  =  1,  NCT!2) 

L  =  NDL!2,JD) 

DO  50,  K  =  1,  3 

ND!K)  =  LND!L,K) 

IF!ND!K).Eq.  2)  KND  =  K 
50  CONTINUE 
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DO  60,  K  =  1,  3 
NL  =  ND(K) 

C  BOUNDARY  CONDITION 

IF(NL.EQ.  1)  THEN 

P(J,1)  =  -UBCOND*FROW(L,KND,K)  +  P(J,1) 

C  UPPER  ROW 

ELSEIF(NL.EQ.  2)  THEN 

B(J,1)  =  FROW(L,KND,K)  +  B(J,1) 

r  TniiTPR  pnu 

ELSEIF((NL.GT.  3).AND.  (NL.LT.  7))  THEN 

C(J,NL-3)  =  FROW(L.KND,K)  +  C(J,NL-3) 

C  ERROR 

ELSE 

WRITE(*,*)  NL, 'ERROR  -  INDEX  EQUALS  3  OR  7 
ENDIF 
CONTINUE 
CONTINUE 

IF((CHAR2.EQ.  'I’).0R.  (CHAR2.EQ.  'i'))  THEN 
WRITECl,*)  I 
DO  72,  M  =  1,  NABC(I,2) 

WRITECl, 1020)  (REAL(B(M,N)),  N  =  1,  NABC(I,2)) 
CONTINUE 

WRITECl,*)  '  ' 

DO  73,  M  =  1,  NABCCI,2) 

WRITECl, 1020)  CREALCCCM,N)),  N=  1,  NABCCI,3)) 
CONTINUE 

WRITECl,*)  ’  ’ 

DO  74,  M  =  1,  NABCCI,2) 

WRITEC 1,1020)  CRIAL(P(M,N)),  N  »  1,  PERND) 

CONTINUE 

WRITECl,*)  '  ’ 

WRITECl,*)  '  ' 

WRITECl,*)  '  ' 

ENDIF 

IFC.  NOT.  CC CHARS.  EQ.  'M' ).  OR.  CCHAR5.  EQ.  ’m’)))  THEN 
CALL  MARCHC I , IMX , NBMX , NABCC 1,1), NABCC 1,2), NABCC 1,3), INORN , SVTC ) 
CALL  ZERO 
ENDIF 
C  BEGIN  LOOPING 

DO  900,  I  =  1,  BIND 
C 

DO  400,  J  =  1,  NDBOT-2 
NOD  *  J  +  NDTOP  +  1 
DO  300,  JD  =  1,  NCTCNOD) 

L  »  NDLCNOD,JD) 

DO  100,  K  =  1,  3 

NDCK)  «  LNDCL,K) 

IFCNDCK).EQ.NOD)  KND  =  K 
100  CONTINUE 

DO  200,  K  =  1,  3 
NL  =  NDCK) 

C  BOUNDARY  CONDITION 

IFCCI.EQ.  1).  AND.  CNL.EQ.  D)  THEN 
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P(J,NL)  =  -UBCOND*FROW(L,KND,K)  +  P(J,NL) 

ELSEIFCNL.EQ.  1)  THEN 

P( J , PERND - T+2 ) =-UBCOND*FROW( L ,KND , K)+P( J , PERND - 1+2 ) 

C  SPECIAL  CASE  FOR  NODE  #2  AND  I  =  1 
ELSE  IF(  ( I .  EQ.  1) .  AND.  (  NL.  EQ.  2  )  )THEN 
A(J,1)  =  FROW(L,KND,K)  +  A(J,1) 

ELSEIFCNL.EQ.  NDTOP)  THEN 

P(J,I)  =  -UBCOND*FROW(L,KND,K)  +  P(J,I) 

ELSEIFCNL.EQ.  NDTOP+1)  THEN 

PC J , PERND- I+l )=-UBCOND*FROWC L ,KND , K)+PC J , PERND- I+l ) 
ELSEIFCNL.EQ.  NDTOT)  THEN 

PC J, I+l)  =  -UBCOND*FROWCL,KND,K)  +  PC J, I+l) 

C  UPPER  ROW 

ELSEIFCNL.lt.  NDTOP)  THEN 

ACJ,NL-1)  =  FROWCL.KND.K)  +  ACJ,NL-1) 

C  LOWER  ROW 

ELSEIFCNL.lt. NDTOT)  THEN 

BCJ,NL-NDT0P-1)  =  FROWCL,KND,K)  +  BC J,NL-NDTOP-l) 

C  ERROR 

ELSE 

WRITEC*,*)  NL, 'ERROR  -  DUE  TO  INDEX  GREATER  THAN  NODE  TOTAL' 

ENDIF 

C 

200  CONTINUE 

300  CONTINUE 

400  CONTINUE 

INDEX  FOR  NEXT  ROW  AND  COMPLETE  B,  C,  P  FILLS 

JJ  =  I  +  1 

WRITEC6,1010)  JJ.TCALL 
CALL  SORTERCJJ.LEL.LAND) 

CALL  FINDERCJJ, LAND, OFFSET) 

IFC.  NOT.  CC CHARS.  EQ.  'M' ).  OR.  CCHAR5.  EQ. 'm')))  THEN 
CALL  BNDC C  J J , BCOND , EO , SURBC , CHAR 1 , ALPHA , BETA , LINE , MODE , 

C  X0RIGIN,Y0RIGIN,CHAR4) 

ENDIF 

CALL  F I LLCJJ,LEL,FROW, ALPHA, BETA) 

C  ESTABLISH  THE  NUMBER  OF  NODES  IN  THE  I+lth  AND  I+2th  ROWS 
IFCJJ.NE.  CBIND+D)  THEN 

NDTOP  =  N0DESCPERND+2-JJ)  +  NODESCJJ)  -  1 
NDBOT  =  NODESCPERND+l-JJ)  +  NODESCJJ+1)  -  1 
NDTOT  =  NDBOT  +  NDTOP 

ELSE 

NDTOP  =  5 
NDBOT  =  2 
NDTOT  =  7 

ENDIF 

C 

DO  800,  J  =  1,  NDTOP-2 
NOD  =  J  +  1 

DO  700,  JD  =  1,  NCTCNOD) 

L  =  NDLCNOD.JD) 

DO  500,  K  =  1,  3 

NDCK)  =  LNDCL.K) 
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IF(ND(K).EQ.NOD)  KND  «  K 
500  CONTINUE 

DO  600,  K  =  1,  3 
NL  =  ND(K) 

C  BOUNDARY  CONDITION 
IFCNL.EQ.  1)  TW^N 

P( J , PERND- +2 )=-UBCOND*FROW( L . KND , K) +P( J , PERND- J J+2 ) 
ELSEIFCNL.  FQ.NDTOP)  THEN 

P(J,JJ)  -UBCOND*FROW(L,KND,K)  +  P(J,JJ) 

ELSEIFC  NL.  I.Q.  NDTOP+1 )  THEN 

P( J , PERND- JJ+1 )=-UBCOND*FROW( L, KND ,K)+P( J , PERND- JJ+1) 

ELSEIFC (JJ.  EQ.  (BIND+1) ).  AND.  (NL. EQ.  NDTOT))  THEN 
C(J,1)  =  FROW(L,KND,K)  +  C(J,1) 

ELSEIFCNL.EQ.  NDTOT)  THEN 

P(J,JJ+1)  =  -UBCOND*FROW(L,KND,K)  +  P(J,JJ+1) 

C  UPPER  ROW 

ELSEIFCNL.  LT.NDTOP)  THEN 

B(J.NL-l)  =  FROW(L,KND,K)  +  B(J,NL-1) 

C  LOWER  ROW 

ELSEIFCNL.  LT.  NDTOT)  THEN 

C(J,NL-NDT0P-1)  =  FROW(L,KND,K)  +  C( J,NL-NDT0P-1) 

C  ERROR 

ELSE 

WRITE (*,*)  NL, 'ERROR  -  DUE  TO  INDEX  GREATER  THAN  NODE  TOTAL' 
ENDIF 
C 

600  CONTINUE 

700  CONTINUE 

800  CONTINUE 

C 

IF( ( CHAR2.  EQ. ' I ' ) . OR, ( CHAR2. EQ. ' i ' ) )  THEN 
WRITECl,*)  JJ 
DO  91,  M  =  1,  NABC(JJ,2) 

WRITECl, 1020)  (REAL(A(M,N)),  N  =  1,  NABC(JJ,1)) 

91  CONTINUE 

WRITECl,*)  '  ' 

DO  92,  M  =  1,  NABC(JJ,2) 

WRITECl, 1020)  (REAL(B(M,N)),  N  =  1,  NABC(JJ,2)) 

92  CONTINUE 

WRITECl,*)  '  ' 

DO  93,  M  =  1,  NABC(JJ,2) 

WRITECl, 1020)  (REAL(C(M,N)),  N  =  1,  NABC(JJ,3)) 

93  CONTINUE 

WRITECl,*)  '  ' 

DO  94,  M  =  1,  NABC(JJ,2) 

WRITECl, 1020)  (REAL(P(M,N)),  N  =  1,  PERND) 

94  CONTINUE 
C 

WRITECl,*)  '  ’ 

WRITECl,*)  '  ' 

WRITECl,*)  '  ' 

ENDIF 

C 

IF(.  NOT.  ((CHARS.  EQ.  'M').0R.  (CHARS. EQ. 'm')))  THEN 
CALL  MARCHC  J J , IMX , NBMX , NABC ( J J , 1 ) , NABC ( J J , 2 ) , NABCC  J J , 3 ) , INORM , 
CSVEC) 
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CALL  ZERO 
END  IF 
C 

900  CONTINUE 
C 

C  LOAD  A,  B  &  P  FOR  I  -  BIND  ♦  2  (BOTTOH) 

C 

J  ■  1 

DO  30,  JD  •  1.  NCT(7) 

L  ■  NDL(7,JD) 

DO  10.  K  »  1.  3 

ND(K)  -  LNTXL.K) 

IF(ND(K).EQ.  7)  KND  •  K 
10  CONTINUE 

DO  20.  K  -  1.  3 
NL  -  KD(K) 

C  BOUNDARY  CONDITION 

IF(NL.EQ.  6)  THEN 

P(J.JJ+1)  ■  -UBCOND*FROW(L.KND.K)  ♦  P(J.JJ*1) 

C  UPPER  ROW 

ELSEIF((NL.GT.  1).  AND.  (NL.  LT.  5))  THEN 

A(J,NL-1)  -  FROW(L,KND.K)  ♦  A(J,SL-l) 

C  LOWER  ROW 

ELSEIF(NL.  EQ.  7)  THEN 

B(J,1)  -  FROW(L.KND,K)  ♦  B(J.l) 

C  ERROR 

ELSE 

WRITEC*,*)  NL, 'ERROR  -  INDEX  EQUAL  TO  1  OR  5* 

ENDIF 

C 

20  CONTINUE 

30  CONTINUE 

C 

IF((CHAR2.EQ. ' I ' ). OR. (CHAR2. EQ. 'i'))  THEN 
WRITECl,*)  TCALL+1 
DO  61,  M  «  1,  NABC( TCALL+1, 2) 

WRITECl, 1020)  (REAL(A(H,N)),  N  «  1,  NABCC TCALL+1 . 1) ) 

61  CONTINUE 

WRITECl,*)  '  ' 

DO  62,  M  «  1,  NABCC TCALL+1, 2) 

WRITECl, 1020)  CREALCBCM,N)),  N  =  1,  NABC C TCALL+1 ,2) ) 

62  CONTINUE 

WRITECl,*)  '  ' 

DO  64,  M  =  1,  NABCC TCALL+1, 2) 

WRITECl, 1020)  CREALCPCM.N)),  N  «  1,  PERND) 

64  CONTINUE 

C 

WRITECl,*)  '  ' 

WRITECl,*)  '  ' 

WRITECl,*)  ’  ' 

ENDIF 

C 

WRITEC6,*)  ’  FINAL  INVERSION' 

IFC.NOT. CCCHAR5.EQ.  'M' ).  OR.  CCHAR5.  EQ.  'm')))  THEN 
CALL  MARCHC  JJ+1 , IMX , NBMX , NABCC  JJ+1 , 1 ) , NABCC  JJ+1 , 2 ) , NABCC  JJ+ 1 , 3 ) , 
CINORM.SVEC) 
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ENDIF 

C 

WRITE(2,*)  'END  END  ’ 

WRITE(6,*)  '  ’ 

WRITE(6,*)  'MINIMUM  AREA  =  ’ .MINAREA 

WRITE(6,*)  'AT  ROW  ' .MINROW,  AND  ELEMENT  NUMBER  '.MINED 

WRITE(6,*)  'MAXIMUM  AREA  =  '.MAXAREA 

WRITE(6,*)  'AT  ROW  ' .MAXROW/  AND  ELEMENT  NUMBER  ' ,MAXEL 
RATIO  =  MAXAREA/MINAREA 
WRITE(6,*)  'AREA  RATIO  =  '.RATIO 
IF(RATIO.GT.  2.  5)  THEN 

WRITE(6,*)  'YOU  SHOULD  CONSIDER  ABORTING  THIS  RUN  AND  ' 
WRITE(6,*)  'LOOKING  AT  THE  MESH  IN  CURVE  DIGITIZER  ' 
WRITE(6,*)  'A  BETTER  METHOD  MAY  BE  AVAILABLE  ' 

ENDIF 

IF((CHAR3.EQ.  'D').OR.  (CHARS.  EQ.  'd'))  THEN 
WRITE(20,1080) 

WRITE(20,1090) 

WRITE(20,1100) 

ENDIF 

C 

1000  FORMAT(  16,’  OUT  OF', 16,’  CALLS  ’) 

1010  FORMATC  16,'  OUT  OF', 16,'  CALLS  ’) 

1020  FORMATC  20(E14.  8 , 1X,E14.  8 , IX) ) 

1030  FORMATC  6X , ' CALL  COMPRS ' ) 

1040  FORMATC  6X, ' CALL  NOBRDR' ) 

1050  FORMATC  6X, 'CALL  PAGE(8. 0,10. 0) ' ) 

1060  FORMATC  6X, ' CALL  AREA2D(5.  0, 7.  0) ' ) 

1070  FORMATC  6X,' CALL  FRAME') 

1080  FORMATC  6X,' CALL  DONEPL') 

1090  FORMATC  6X,'STOP') 

1100  FORMATC  6X,'END') 

C 

RETURN 

END 

C 


C 

SUBROUTINE  FILL( I , LEL , FROW , ALPHA , BETA) 

C 

COMPLEX  F(3,3),  FROWC 100,3,3) ,  ALPHA,  BETA,  A,  B 
REAL  AREA,  MINAREA,  MAXAREA 
REAL  NDP(200,2) 

INTEGER  I,  J,  K,  L,  LEL,  LND(0: 200,3) ,  NDL(200,4),  NCT(200) 

INTEGER  MINROW,  MINED,  MAXROW,  MAXEL,  M 

CHARACTER*!  CHAR2,  CHAR3,  CHAR4 

C0MM0N/BLK4/LND , NDL , NCT 

C0MM0N/BLK5/NDP 

COMMON/BLK7/MINAREA , MINROW .MINED , MAXAREA .MAXROW , MAXEL , AREA 
COMMON/BLK9/CHAR2,  CHAR3,  CHAR4 
C 

IF((CHAR2.EQ.  '  I '  ).  OR.  (CHAR2.  EQ.  'l'))  THEN 
WRITECll,*)  'ROW  NUMBER  =  ',! 

ENDIF 

DO  30,  J  =  1,  LEL 

IF((CHAR4.EQ. ' U' ) . OR. (CHAR4.  EQ. 'u'))  THEN 
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A  *  ALPHA 
B  »  BETA 

ELSE 

IF((J.  EQ.  1).  OR.  (J.EQ.  2).  OR.  (J.  EQ.  LEL-l).  OR.  ( J.  EQ.  LEL)) 
C  THEN 

A  -  (1.0, 0.0) 

B  =  (1.0, 0.0) 

ELSE 

A  «  ALPHA 
B  ■  BETA 

ENDIF 

ENDIF 

CALL  VARINT(J,F,A.B,AREA,LND) 

IF((J.GT.  2).AND.(J.LT.LEL-1))  THEN 
IF(AREA.  LT.MINAREA)  THEN 
MINAREA  «  AREA 
MINROW  «  I 
MINEL  «  J 

ELSEIF(AREA.GT.MAXAREA)  THEN 
MAXAREA  «  AREA 
MAXROW  =  I 
MAXEL  =  J 

ENDIF 

ENDIF 

DO  20,  K  =  1,  3 

WRITE(2,*)  NDP(LND(J,K),1),  NDP(LND( J,K) ,2) 

DO  10,  L  =  1,  3 

FROW(J,K,L)  =  F(K,L) 

10  CONTINUE 

IF((CHAR2.EQ.  'I’).0R.(CHAR2.EQ.  ’i’))  THEN 

WRITE(11,1000)  J,K,(REAL(F(K,L)),  L«  1,3) 

ENDIF 

20  CONTINUE 

IF((CHAR2.EQ.  ' I ' ).  OR. (CHAR2. EQ.  ’i’))  THEN 
WRITEdl,*)  '  ' 

ENDIF 

WRITE(2,*)  NDP(LND(J,1),1),  NDP(LND( J, 1) ,2) 

WRITE(2,*)  '  999990  999990  ' 

C  DISSPLA  PROGRAM  GENERATION 

IF((CHAR3.EQ.  'D' )•  OR. (CHAR3.  EQ.  ’d’))  THEN 

WRITE(20,1010)  NDP(LND(J,1),1),  NDP(LND( J, 1) ,2) 
WRITE(20,1020)  NDP(LND(J,2) , 1) ,  NDP(LND( J,2) ,2) 
WRITE(20,1020)  NDP(LND( J,3) ,1) ,  NDP(LND( J,3) ,2) 
WRITE(20,1020)  NDP(LND( J, 1) , 1) ,  NDP(LND( J, 1) ,2) 

ENDIF 

C 

30  CONTINUE 

C 

1000  FORMAT( IX , 2( 13 , 2X) , 3X , 3(F8. 5 , 2X) ) 

1010  F0RMAT(6X.'CALL  STRTPT( ' ,F8.  5 , ' , ' ,F8.  5 , ’ ) ' ) 

1020  F0RMAT(6X,'CALL  C0NNPT( ' ,F8.  5, ' , ' ,F8.  5, ' )’ ) 

C 

RETURN 

END 

C 
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c 


C 

C 


40 


60 


SUBROUTINE  ZERO 

ZERO  A,  B,  C,  AND  P  MATRICES 

COMPLEX  A(50,50),  B(50,50),  C(50,50),  P(50,100) 
INTEGER  J,  K,  L 

C0MM0N/BLK8/A,B,C,P 


DO  50,  J  =  1,  50 
DO  40,  K  =  1, 
A(J,K)  = 
B(J,K)  = 
C(J,K)  = 
CONTINUE 
DO  60,  L  =  1, 
PCJ.L)  = 
CONTINUE 


50 

CMPLX(0.0,0.  0) 
CMPLX(0.  0.0.  0) 
CMPLX(0.  0,0.0) 

100 

CMPLX(0.  0,0.0) 


50  CONTINUE 
RETURN 
END 


SUBROUTINE  BNDC ( I , BCOND , EO , SURBC , CHARI , ALPHA , BETA , LINE , MODE , 
CX0RIGIN,Y0RIGIN,CHAR4) 

BOUNDARY  CONDITION  FILL  FOR  A  PLANE  WAVE 

EO  =  FIELD  STRENGTH 

KO  =  WAVE  NUMBER  (2*PI /WAVELENGTH) 

BOUNDARY  CONDITION  FILL  FOR  A  CYLINDRICAL  BOUNDARY  CONDITION 
EO  =  FIELD  STRENGTH 

MODE  =  MODE  NUMBER  FOR  CYLINDRICAL  BOUNDARY  CONDITIONS 

COMPLEX  SURBC(IOO),  VALUE(50),  BCOND(IOO),  ALPHA,  BETA,  LINE(50) 
COMPLEX  K  RAK 

MAL  NDP(200,2),  EO,  KR,  KI,  XORIGIN,  YORIGIN,  RA,  RB 

INTEGER  I,  J,  NDTOP,  NDBOT,  NDTOT,  NODES(200),  PERND,  BIND,  MODE 

CHARACTER*!  CHARI,  CHAR4 

C0MM0N/BLK2/PERND , BIND 

C0MM0N/BLK3/N0DES 

C0MM0N/BLK5/NDP 

IF((CHAR1.EQ.  'P').OR.  (CHARI.  EQ.  ’p’))  THEN 

IF((CHAR4.EQ.  ' U *  ) .  OR.  ( CHAR4.  EQ.  'u'))  THEN 
KR  =  REAL(CSQRT(BETA)) 

KI  =  ABS(AIMAG(CSQRT(BETA))) 

ELSE 

KR  =  1.  0 
KI  =  0.  0 

ENDIF 

CHECK  ON  WHETHER  WE  ARE  USING  ER  =  ALPHA  OR  BETA 
IFCI.EQ. 1)  THEN 

BCOND(l)  =  E0*EXP(KI*NDP(1,2))*CMPLX(C0S(KR*NDP(1,2)), 
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C  SIN(KR*NDP(1,2))) 

BCOND(PERND)  =  EO*EXP(KI*NDP(3,2))*CMPLX(COS(KR*NDP(3,2)) , 

C  SIN(KR*NDP(3,2))) 

BCOND(2)  =  E0*EXP(KI*NDP(7,2))*CMPLX(C0S(KR*NDP(7,2)), 

C  SIN(KR*NDP(7,2))) 

VALUE(l)  =  E0*EXP(KI*NDP(2,2))*CMPLX(C0S(KR*NDP(2,2)), 

C  SIN(KR*NDP(2,2))) 

WRITEdO,*)  BCOND(l) 

WRITE(10,1000)  VALUE(l) 

WRITEdO,*)  ’  ' 

SURBCd)  =  VALUECl) 

ELSEIFd.LE.  BIND)  THEN 

NDTOP  =  NODES(I)  +  NODES(PERND-I+2)  -  1 
NDBOT  =  NODESd+1)  +  NODES(PERND-I+l)  -  1 
NDTOT  =  NDTOP  +  NDBOT 

BCOND(PERND-I+l)  =  EO*EXP(KI*NDP(NDTOP+1,2))*CMPLX(COS(KR* 

C  NDP( NDTOP+1 , 2 ) ) , SIN( KR*NDP( NDTOP+1 , 2 ) ) ) 

BCOND(I+l)  =  EO*EXP(KI*NDP(NDTOT,2))*CMPLX(COS(KR*NDP(NDTOT, 
C  2)),SIN(KR*NDP(NDTOT,2))) 

WRITEdO,*)  BCOND(PERND-I+2) 

DO  10,  J  =  2,  NDTOP-1 

VALUE(J)  =  EO*EXP(KI*NDP(J,2))*CMPLX(COS(KR*NDP(J,2)), 

C  SIN(KR*NDP(J,2))) 

IF(J.EQ.  2)  THEN 

SURBC(PERND-I+2)  =  VALUE(2) 

ELSEIFCJ.EQ.  (NDTOP-D)  THEN 
SURBC(I)  =  VALUE(NDTOP-l) 

ENDIF 

10  CONTINUE 

LINE(I)  =  VALUE(( NDTOP+1 )/2) 

WRITEdO,  1000)  (VALUE(J),  J  =  2,  NDTOP-1) 

WRITEdO,*)  BCOND(I) 

WRITEdO,*)  '  ' 

C 

ELSEIFCI.EQ.  BIND+1)  THEN 

BCONDCI+l)  =  E0*EXP(KI*NDP(6,2))*CMPLX(C0S(KR*NDP(6,2)), 

C  SIN(KR*NDP(6,2))) 

VALUE(2)  =  E0*EXP(KI*NDP(2,2))*CMPLX(C0S(KR*NDP(2,2)), 

C  SIN(KR*NDP(2,2))) 

VALUE(3)  =  EO*EXP(KI*NDP(3,2))*CMPLX(COS(KR*NDP(3,2)), 

C  SIN(KR*NDP(3,2))) 

VALUE(4)  =  E0*EXP(KI*NDP(4,2))*CMPLX(C0S(KR*NDP(4,2)), 

C  SIN(KR*NDP(4,2))) 

WRITEdO,*)  BC0ND(I+2) 

WRITEd0,1000)  VALUE(2),  VALUE(3),  VALUE(4) 

SURBC(BIND+3)  =  VALUE(2) 

SURBC(BIND+1)  =  VALUE(4) 

LINECI)  =  VALUE(3) 

WRITEdO,*)  BCONDd) 

WRITEdO,*)  '  ' 

VALUE(2)  =  EO*EXP(KI*NDP(7,2))*CMPLX(COS(KR*NDP(7,2)), 

C  SIN(KR*NDP(7.2))) 

WRITEdO,  1000)  VALUE(2) 

WRITEdO,*)  BCONDCI+l) 

SURBC(BIND+2)  =  VALUE(2) 
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ENDIF 

C 

ELSEIFC (CHARI.  EQ.  'C').OR.  (CHARI.  EQ.  'c'))  THEN 
C 

IFd.EQ.  1)  THEN 

RA  =  SQRT((NDP(2,1)-X0RIGIN)**2+(NDP(2,2)“Y0RIGIN)**2) 

RB  =  SQRT((NDP(1,1)-X0R1GIN)**2+(NDP(1,2)-Y0RIGIN)**2) 

WRITE (*,*)  BETA 
K  =  CSQRT(BETA) 

RAK  =  RA*K 

CALL  CYLBC( RA , RB , RAK , NDP( 1,1). NDP( 1 , 2 ) , XORIGIN , YORIGIN , EO , 

C  BCOND(l),SURBC(l),K) 

CALL  CYLBC(RA,RB, RAK. NDP(3,1),NDP(3, 2), XORIGIN, YORIGIN, EO, 

C  BCOND( PERND) , SURBC( PERND) ,K) 

CALL  CYLBC( RA , RB , RAK , NDP( 7 , 1 ) , NDP( 7 , 2 ) , XORIGIN , YORIGIN , EO , 

C  BCOND(2),SURBC(2),K) 

WRITEdO,*)  BCOND(l),  SURBC(l) 

WRITEdO,*)  '  ' 

ELSEIFd.LE.  BIND)  THEN 

NDTOP  =  NODESd)  +  NODES( PERND- 1+2)  -  1 
NDBOT  =  NODES(I+l)  +  NODES( PERND- I+l)  -  1 
NDTOT  =  NDTOP  +  NDBOT 

CALL  CYLBC(RA,RB,RAK,NDP(NDTOP+l,l) ,NDP(NDTOP+l,2) .XORIGIN, 
C  YORIGIN, EO,BCOND(PERND-I+l),SURBC(PERND-I+l),K) 

CALL  CYLBC( RA , RB , RAK , NDP( NDTOT , 1 ) , NDP( NDTOT , 2 ), XORIGIN , 

C  YORIGIN, EO ,BCOND( I+l) ,SURBC( I+l) ,K) 

WRITEdO,*)  BCOND( PERND -1+2),  SURBC(PERND-I+2) 

WRITEdO,*)  BCOND(I),  SURBC(I) 

WRITEdO,*)  '  ' 

ELSEIFd.EQ.  BIND+1)  THEN 

WRITEdO,*)  BCOND(I+2),  SURBC(I+2) 

WRITEdO,*)  BCOND(I),  SURBC(I) 

WRITEdO,*)  '  ' 

CALL  CYLBC(RA,RB,RAK,NDP(6,1),NDP(6,2),X0RIGIN, 

C  YORIGIN, EO,BCOND( I+l), SURBC( I+l), K) 

WRITEdO,*)  BCONDd+l),  SURBC(I+1) 

WRITEdO,*)  '  ' 

ENDIF 

C 

ELSE 

RETURN 

ENDIF 

C 

1000  F0RMAT(1X,50(E14.  8,2X,E14. 8,2X)) 

C 

RETURN 

END 


SUBROUTINE  CYLBC(RA,RB, RAK, X.Y, XORIGIN, YORIGIN, EO,BC,PSI,K) 

C 

COMPLEX  SQRTMl,  JORB,  JIRB,  JORAK,  JIRAK,  HORA,  HIRA,  HORB,  HIRB 
COMPLEX  DELTAN,  AN,  PSI,  JORA,  JIRA,  K,  RAK,  BC 
REAL  PI,  XORIGIN,  YORIGIN,  X,  Y,  RA,  RB,  EO,  PHI 
C 

PI  =  4.  0*ATANd.O) 
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SQRTMl  =  CMPLX(0.  0,1. 0) 

PHI  =  ATAN2(X  -  XORIGIN,  Y  -  YORIGIN) 

CALL  BES1(CMPLX(RA.0.  0),J0RA,J1RA) 

CALL  BES1(CMPLX(RB,0.  0) , JORB , JIRB) 

CALL  BES1(RAK,J0RAK,J1RAK) 

CALL  HAN1(CMPLX(RA,0.  0),H0RA,H1RA) 

CALL  HAN1(CMPLX(RB,0. 0),H0RB,H1RB) 

C 

DELTAN  =  J1RB*(J1RAK*(H0RA-H1RA/RA)-K*(J0RAK-J1RAK/RAK)*H1RA)- 
CH1RB*(J1RAK*(J0RA-J1RA/RA)-K*(J0RAK-J1RAK/RAK)*J1RA) 

AN  =  -2.  0*SQRTM1/(PI*RA*DELTAN) 

BC  =  EO*COS(PHI) 

PSI  =  AN*J1RAK*BC 
C 

RETURN 

END 


SUBROUTINE  BES1(Z,J0,J1) 

Computing  Bessel  Functions  for  n  =  0,  1  with 
Complex  Argument  Z.  Direct  Power  Series  Method  for 
CABS(Z)  .  LE.  6  and  Hankel's  Asymptotic  Formula  for 
CABS(Z)  .GT.  6. 

Written  11/5/87  by  M.  A.  Morgan 
INTEGER  M,M2 

REAL  C( 34) , DM , F( 34 ) , GO , P( 34) , Pi , P2 

COMPLEX  Z,Z2,Z3,Z4,J0,J1,AM,CL,P0,P1,Q0,Q1,C0,C1,S0,S1 
Pi=3. 1415927 
P2=2. 0/PI 

IF(CABS(Z).LE.  6.  0)  THEN 

Utilizing  the  Direct  Power  Series  Method 

G0=  1.781072 
Z2=0. 5*Z 
CL=CL0G(G0*Z2) 

Computing  F(m)  =  m  !  and  P(m)  =  1  +  1/2  +  1/3  +  ....+  1/m 

F(l)=1.0 
P(l)=l.  0 
DO  11  M=2,34 

F(M)=M*F(M-1) 

P(M)=P(M-1)+1.  0/M 
1  CONTINUE 

Computing  Power  Series  Coefficients 
C 

DM=-1.  0 
DO  22  M=l,34 

C(M)=DM/(F(M)*F(M)) 

DM=-DM 

22  CONTINUE 

C 
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Computing  JO  and  J1 

J0a(l.  ,0.  ) 

J1»(0. ,0. ) 

M=0 

3  M*M+1 

M2=2*M 

AM*C(M)*(Z2**M2) 

J0=J0+AM 

IF((CABS(AM).GT.  1.0E-10).AND.  (M.  LT.  34))  GO  TO  33 

J1=J1/Z2 

return 

ELSE 

Hankel'  Asymptotic  Formula  (Abram.  &  Stegun  p.  364) 

Z2=Z*Z 

Z3=Z*Z2 

Z4=Z*Z3 

P0=1. 0-.  0703125/Z2+.  1121521/Z4 
Q0=-. 125/Z+.  0732422/Z3 
Pl=l. 0+. 1171875/Z2-.  1441956/Z4 
Ql=.  375/Z-.  10253906/Z3 
CO=CCOS(Z-.  25*PI) 

S0=CSIN(Z-.  25*PI) 

C1=CC0S(Z-.  75*PI) 

S1=CSIN(Z-.  75*PI) 

AM=CSQRT(P2/Z) 

J0=AM*(P0*C0-Q0*S0) 

J1=AM*(P1*C1-Q1*S1) 

ENDIF 

RETURN 

END 

SUBROUTINE  HAN1(Z,H0,H1) 


Computing  Hankel  Functions  for  n  =  0,  1  with 
Complex  Argument,  Z.  Direct  Power  Series  Method  for 
CABS(Z)  . LE.  5  and  Hankel's  Asymptotic  Formula  for 
CABS(Z)  .GT.  5.  Written  11/6/87  by  M. A.  Morgan 

INTEGER  M,M2 

REAL  C( 34) , DM , F( 34 ) , GO , P ( 34 ) , Pi , P2 

COMPLEX  Z,Z2,Z3,Z4,J0,Jl,Y0,yi,AM,CL,P0,Pl,Q0,Ql 

COMPLEX  E0,El,X0,Xl,H0,Hl,j 

PI=3. 1415927 

P2=2.  0/PI 

jKO.  ,1.  ) 

IF(CABS(Z).LE.  5.  0)  THEN 

Direct  Power  Series  Method 

G0=  1.78072 
Z2=0. 5*Z 


no 


CL=CL0G(G0*Z2) 

C 

C  Computing  F(m)  =  m  !  and  P(m)  =  1  +  1/2  +  1/3  +  - +  1/m 

C 

F(l)=1.0 
P(l)=l.  0 
DO  11  M=2,34 

F(M)=M*F(M-1) 

P(M)=P(M-1)+1.0/M 
11  CONTINUE 

C 

C  Computing  Power  Series  Coefficients 
C 

DM=-1.  0 
DO  22  M=l,34 

C(M)=DM/(F(M)*F(M)) 

DM=-DM 

22  CONTINUE 

C 

C  Computing  JO  and  J1 
C 

J0=(1.  ,0.) 

J1=(0. ,0. ) 

M=0 

33  M=M+1 

M2=2*M 

AM=C(M)*(Z2**M2) 

J0=J0+AM 

J1=J1-M*AM 

IF((CABS(AM).GT.  1.0E-10).AND.(M.  LT.  34))  GOTO  33 
J1=J1/Z2 
C 

C  Computing  YO  and  Yl 
C 

M=0 

Y0=CL*J0 

Yl=Z2*CL*Jl-0.  5*J0 
44  M=M+1 

M2=2*M 

AM=C(M)*P(M)*(Z2**M2) 

Y0=Y0-AM 

Y1=Y1+M*AM 

IF((CABS(AM).GT, 1.0E-10).AND. (M. LT. 34))  GO  TO  44 

Y0=P2*Y0 

Y1=P2*Y1/Z2 

H0=J0-j*Y0 

Hl=Jl-j*Yl 

RETURN 

ELSE 

C 

C  Hankel'  Asymptotic  Formula  (Abram.  &  Stegun  p.  364 
C 

Z2=Z*Z 

Z3=Z*Z2 

Z4=Z*Z3 

P0=1.  0-.  0703125/Z2+. 1121521/Z4 


in 
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Q0=-.  125/Z+. 0732422/Z3 
Pl=1.0+.  1171875/Z2-.  1441956/Z4 
Ql=.  375/Z-. 10253906/Z3 
XO=(Z-. 25*PI) 

X1=(Z-. 75*PI) 

EO=CEXP(-j*XO) 

El=CEXP(-j*Xl) 

AM=CSQRT(P2/Z) 

HO=AM*(PO-J*QO)*EO 

Hl=AM*(Pl-j*Ql)*El 

ENDIF 

RETURN 

END 

SUBROUTI NE  MARCH( I , I MX , NBMX . NA , NB . NC , I NORM , S VTC ) 

THIS  ROUTINE  PERFORMS  THE  RICCATl  TRANSFORM  FIRST 
SWEEP,  GENERATING  AND  STORING  ON  DISK  1  RMAT 
AND  SVEC  (FOR  EACH  MODE)  AT  EACH  FORWARD  STEP. 

COMPLEX  RMAT(50.S0) .SVEC(50,100) 

COMPLEX  A(SO,50).B(50.50),C(50,50),P(50,100) 

COMPLEX  D(50).SUM,DET 
REAL  COND.DMAG 

INTEGER  I.J.K.L.NA.NB.NC.PERND.BIND.IMX.INORM 
INTEGER  NBMX 

CHARACTER*!  CHAR2,  CHAR3,  CHAR4 

COMMON/ BLK2 / PERND .BIND 
COMMON/ BLK8/ A, B,C,P 
COMMON/BLK9/CHAR2 ,  CHAR3,  CHAR4 

LOADING  THEN  INVERTING  (B+A*RMAT) 

USING  MINIMUM  MEMORY  SINGLE  MATRIX  TECHNIQUE 
DATE  1/29/80  FOR  THIS  CHANGE 

SKIPPING  FIRST  A*R  (WHEN  A  «  0,  ZERO  R  MATRIX) 
IFCI.EQ.  1)  THEN 

DO  10,  J  «  1,  NB 

DO  10,  K  »  1,  NB 

RMAT(J,K)  ■  (0.  ,0.  ) 

0  CONTINUE 

RMAT  =  A*R 
ELSE 
C 

IF((CHAR2.  EQ.  ’l’).OR.  (CHAR2.EQ.  ’i’))  THEN 
WRITE(19,*)  '  OLD  R  MATRIX' 

DO  11,  J  »  1,  5 

WRITE(19,1000)  (REAL(RMAT(J,K)),  K  -  1,  5) 

11  CONTINUE 
WRITE(19,*)  '  ' 

WRITE(19,*)  '  A  MATRIX’ 

DO  12,  J  =  1,  5 

WRITE(19,1000)  (REAL(A(J,K)),  K  »  1,  5) 

12  CONTINUE 
WRITE(19,*)  '  ' 

ENDIF 

C 


II2 


noon 


20 


30 

C 


13 


C 

C 


14 


C 

C 


40 

C 


15 


C 

C 


DO  30,  K  =  1,  NB 

DO  20.  J  =  1,  NB 

D(J)  *  (0.  .0.  ) 

DO  20,  L  =  1,  NA 

D(J)  =  D(J)  +  A(J,L)*RMAT(L,K) 

CONTINUE 

DO  30.  J  »  1,  NB 

RMAT(J,K)  =  D(J) 

CONTINUE 

IF((CHAR2.EQ.  ’I').0R.  (CHAR2.EQ.  'i'))  THEN 
WRITE! 19,*)  ’  NEW  R  MATRIX’ 

DO  13,  J  ■  1,  5 

WRITE! 19, 1000)  !R£AL!RMAT!J,K)).  K  =  1,  5) 
CONTINLX 
WRITE!19,*)  ’  ' 

ENDIF 

ENDIF 

IF!!CHAR2.  EQ. ' l' ). OR.  !CHAR2. EQ.  ' l' ))  THEN 
WRITE! 19,*)  '  B  MATRIX* 

DO  14,  J  =  1,  5 

WRITE! 19, 1000)  !REAL!B!J,K)),  K  =  1,  5) 
CONTINUE 
WRITE! 19,*)  '  ' 

ENDIF 

RMAT  *  B  +  RMAT 

DO  40.  J  »  1,  NB 

DO  40,  K  =  1,  NB 

RMAT!J,K)  »  RMAT!J,K)  +  B!J,K) 

CONTINUE 

IF!!CHAR2. EQ.  T ' ).  OR.  !CHAR2.  EQ. ’i’))  THEN 
W’RITE!19,*)  '  NTWEST  R  MATRIX’ 

DO  15,  J  =  1,  NB 

WRITE!9,1000)  ! REAL! RMAT! J.K)),  K  =  1,  NB) 

WRITE! 19, 1000)  !REAL!RMAT!J,K)),  K  =  1,  NB) 
CONTINUE 
WRITE!9,*)  ’  ' 

WRITE! 19,*)  ’  ’ 

ENDIF 

INVERTING  THE  MATRIX  !B  +  A*R) 

CALL  CSMINV! RMAT , NBMX , NB , DET , COND , INORM) 

DMAG  =  CABS! DET) 

WRITE!*,*)  I, NBMX, NB, DMAG, COND 


COMPUTING  THE  NEW  S -VECTORS 


SKIPPING  FIRST  A*S  !A  =  0) 

IF  !I.EQ.  1)  THEN 
CONTINUE 

ELSE 

DO  70.  K  =  1,  PERND 
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DO  60,  J  »  1,  NB 

D(J)  »  (0.  0,0.  0) 

DO  60,  L  =  1,  NA 

D(J)  *  D(J)  +  A(J,L)*SVEC(L,K) 

60  CONTINUE 

DO  70,  J  «  1,  NB 

P(J,K)  =  P(J,K)  -  D(J) 

70  CONTINUE 

ENDIF 

C  FINAL  SVEC  MULTIPLICATION 

DO  100,  K  =  1,  PERND 
DO  90,  J  =  1,  NB 

D(J)  =  (0.  0,0.  0) 

DO  90,  L  =  1,  NB 

D(J)  =  D(J)  +  RMAT(J,L)*P(L,K) 

90  CONTINUE 

DO  100,  J  =  1,  NB 

SVEC(J,K)  =  D(J) 

100  CONTINUE 

C 

C  STORING  I+l  SVEC  ON  DISK  7 

DO  110  J  =  1,  NB 

WRITE(7)  (SVEC(J,K),  K  =  1,  PERND) 

110  CONTINUE 
C 

C  FINAL  RMAT  MULTIPLICATION 

IF(I.EQ. IMX)  RETURN 
DO  130,  J  -  1,  NB 

DO  120,  K  =  1,  NC 

D(K)  =  (0.  0,0.  0) 

DO  120,  L  =  1,  NB 

D(K)  =  D(K)  -  RMAT(J,L)*C(L,K) 

120  CONTINUE 

DO  130,  K  =  1,  NC 

RMAT(J,K)  =  D(K) 

130  CONTINUE 
C 

C  STORING  I+l  RMAT  ON  DISK  7 
DO  140,  J  =  1,  NB 
WRITE(7)  (RMAT(J,K),  K  =  1,  NC) 

140  CONTINUE 
C 

1000  F0RMAT(10(E9.  3,1X)) 

C 

RETURN 

END 

C 

C 

SUBROUTINE  SWEEP( IMX , NABC , SURBC , LINE , CHARI , U , BCOND , PS I , ANS , CHARS ) 
C  THIS  ROUTINE  PERFORMS  THE  RICCATI  TRANSFORM  BACKSWEEP 

C  FROM  I=IMX  TO  1=1,  RECALLING  RMAT  AND 

C  SVEC  FROM  DISK  7  AT  EACH  BACKSTEP  TO  FORM 

C  THE  NODE  VECTORS  PSI,  THEN  STORING  THESE  ON 

C  DISK  8  FOR  EACH  APPLIED  DIRICHLET  B.  C. 

COMPLEX  RMAT(50, 100) ,PSI( 50, 100) ,SVEC(50 , 100) ,D( 100) ,TEMP 
COMPLEX  SURBC( 100) ,ANS( 100) ,LINE(50) ,ANSB(50) ,U( 100,100) 
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COMPLEX  BCOND(IOO) 

REAL  ERRORP.ERRORD.TERRN.TERRD.ATSERR 
INTEGER  I,J,K,L,NABC( 100,3) . IMX.PERXD, BIND 
CHARACTER*!  CHARI,  CHARS 
C 

COMMON/ BLK2 / PERND .BIND 
C 

IF((CHAR5.EQ.  'M' ).  OR.  (CHARS.  EQ.  ‘m’))  THEN 
RETURN 

ENDIF 

C 

C  INITIAL  DISK  READ  AT  IMX  (R  ■  0,  NOT  WRITTEN,  ■>  S  IS  READ  FIRST) 
WRITE(*,*)  '  ' 

DO  90,  I  «  IMX.  1.  -1 

WRITE(*,10S0)  (IMX-I+l),  IMX 
IFd.EQ.  IMX)  THEN 

DO  10.  J  ■  NABC(I,2),  1,  -1 
BACKSPACE  7 

READ(7)  (PSKJ.K),  K  -  1,  PERND) 

BACKSPACE  7 

10  CONTINUE 

DO  20.  J  »  1.  NABC(I,2) 

DO  IS,  K  «  1,  PERND 

U(IMX.K)  «  PSI(J,K) 

15  CONTINUE 

20  CONTINUE 

C  SUBSEQUENT  DISK  READS 

ELSE 

C  READ  R  MATRIX 

DO  30,  J  =  NABC(I,2),  1,  -1 
BACKSPACE  7 

READ(7)  (RMAT(J,K),  K  «  1,  NABC(I,3)) 

BACKSPACE  7 

30  CONTINUE 

C  READ  S  VECTOR 

DO  40,  J  =  NABC(I,2),  1,  -1 
BACKSPACE  7 

READ(7)  (SVEC(J,K),  K  =  1,  PERND) 

BACKSPACE  7 

40  CONTINUE 

C  MULTIPLY  RMAT  =  RMAT*PSI 

DO  60.  J  =  1,  NABC(I,2) 

DO  SO,  K  =  1,  PERND 

D(K)  =  (0.0, 0.0) 

DO  SO,  L  =  1,  NABC(I,3) 

D(K)  =  D(K)  +  RMAT(J,L)*PSI(L,K) 

50  CONTINUE 

DO  60,  K  =  1,  PERND 
RMAT(J,K)  =  D(K) 

60  CONTINUE 

C  PSI  =  RMAT  +  SVEC 

DO  70,  J  =  1,  NABC(I,2) 

DO  70,  K  =  1,  PERND 

PSI(J,K)  =  RMAT(J,K)  +  SVEC(J,K) 

70  CONTINUE 

ANSB(I)  =  (0.0,0.  0) 
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DO  80,  J  =  1,  NABC(I.2) 

DO  75,  K  =  1,  PERND 

IF((J.EQ.  NABC(I,2)).OR.  (I.EQ.  1))  THEN 
U(I,K)  =  PSI(J,K) 

ELSEIFCJ.EQ. 1)  THEN 

U(2*IMX-I,K)  =  PSI(J,K) 

ELSEIFCJ.EQ. (NABC( I, 2)+l)/2)  THEN 

ANSB(I)  =  ANSB(I)  +  PSI( J,K)*BCOND(K) 

ENDIF 

75  CONTINUE 

80  CONTINUE 

ENDIF 

90  CONTINUE 

WRITE(8,*)  '  ' 

C 

DO  98,  J  =  1,  PERND 

ANS(J)  =  (0.0,0.  0) 

DO  94,  L  =  1,  PERND 

ANS(J)  =  ANS(J)  +  U(J,L)*BCOND(L) 

94  CONTINUE 

98  CONTINUE 

C 

TERRN  =  0.  0 
TERRD  =  0.  0 
DO  100,  I  =  1,  PERND 

ERRORP  =  (CABSCANSd)  -  SURBC(I)))**2 
ERRORD  «  (CABS(SURBC(I)))**2 
TERRN  =  TERRN  +  ERRORP 
TERRD  =  TERRD  +  ERRORD 

WRITE(13,1020)  I,  BCOND(I),  SURBC(I),  ANS(I),  ERRORP 
100  CONTINUE 

WRITE(13,*)  '  ' 

ATSERR  =  (TERRN/TERRD)**0. 5 
WRITEC 13,1030)  ATSERR 

WRITEC*,*)  'RMS  ERROR  (FOR  THE  PERIMETER)  =  ', ATSERR 
WRITEC 13,*)  ’  ' 

C 

IF((CHAR1.EQ.  'C').OR.  (CHARI.  EQ.  'c'))  THEN 
RETURN 

ELSE 

TERRN  =  0.  0 
TERRD  =  0.  0 
DO  110,  I  =  2,  BIND+1 

ERRORP  =  (CABS(ANSB(I)  -  LINE(I)))**2 
ERRORD  =  (CABS(LINE(I)))**2 
TERRN  =  TERRN  +  ERRORP 
TERRD  =  TERRD  +  ERRORD 

WRITEC 13, 1025)  (I-l),  LINE(I),  ANSB(I),  ERRORP 
110  CONTINUE 

WRITEC 13,*)  '  ' 

ATSERR  =  CTERRN/TERRD)**0. 5 
WRITEC 13, 1040)  ATSERR 

WRITEC*,*)  'RMS  ERRROR  (FOR  BISECTION  SEGMENT)  =  '.ATSERR 

ENDIF 

C 

DO  120,  1=1,  PERND 
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WRITE(30,1060)  (U(I,J),  J  =  1,  PERND) 

120  CONTINUE 
C 

1020  F0RMAT(1X,I3,2X,3(E14.  8,1X,E14.  8,3X),F10.  6) 

1025  F0RMAT(1X,I3,2X,4(E14.  8,2X),F10. 6) 

1030  FORMATCIX,’  RMS  ERROR  (FOR  THE  PERIMETER)  =  ' ,F12.  6) 

1040  FORMATdX,'  RMS  ERROR  (FOR  BISECTION  SEGMENT)  =  ’  ,F12.  6) 

1050  FORMATdX, 1 3.'  TOTAL  BACKSWEEP  ROWS  OUT  OF  ',13,'  COMPLETED') 
1060  FORMAT(1X,100(E8.  2,E8.  2,2X)) 

C 

RETURN 

END 

C 

CC 

C 

SUBROUTINE  CSMINV(A,NDIM,N, DETERM , COND , INORM) 

C 

C  INORM  -  FLAG  TO  NORMALIZE  COLUMNS  AND  ROWS  OF  MATRIX  A 
C  MATRIX  NORMALIZATION  BY  M.  A.  MORGAN 

C  APRIL  24,1978 

C 

C  A  -  MATRIX  TO  INVERT 

C  NDIM 

C  N  - 

C  DETERM  -  DETERMINATE  OF  A 

C  COND  -  CONDITION  NUMBER  OF  A 

C  INORM  -  INTEGER  NORMALIZATION  FLAG 

C 
C 

COMPLEX  A(50,50),PIVOT(50), AMAX ,T , SWAP , DETERM , U 

INTEGER  I , J , K , L , IPIV0T( 50 ) , INDEX( 50 , 2 ) , IROW , ICOLUM , LI , JROW 

INTEGER  JCOLUM,N, INORM 

REAL  TEMP , ALPHA( 50), COL( 50 ) , ROW( 50 ) , AJK , SUMAXA , SUMROW , SUMAXI 
C 

IF(NDIM.GT.  50)  THEN 

WRITE(*,*)  '  ERROR  IN  INVERTION  CALL...  DIMENSION  >  50  ' 
STOP 

ENDIF 

IF(N.  GT.NDIM)  THEN 

WRITE(*,*)  '  ERROR  IN  INVERTION  CALL.  . .  N  >  MAX  DIM.  ' 
STOP 

ENDIF 

IF( INORM.  NE.  1)  GO  TO  7 
DO  3  K  =  1,  N 

COL(K)  =  0.0 
DO  1  J  =  1,N 

AJK  =  CABS(A(J,K)) 

IF(AJK. GT.COL(K))  COL(K)  =  AJK 

1  CONTINUE 

DO  2  J  =  1,  N 

2  A(J,K)  =  A(J,K)/C0L(K) 

3  CONTINUE 

C  ROW  NORMALIZING 

DO  6  J  =  1,  N 

ROW(J)  =0.0 
DO  4  K  =  1,  N 


-  INPUT/OUTPUT 

-  INPUT 

-  INPUT 

-  OUTPUT 

-  OUTPUT 

-  INPUT 


AJK  =  CABS(A(J,K)) 

IF(AJK.GT.ROW(J))  ROW(J)  *  AJK 

4  CONTINUE 

DO  5  K  =  1,  N 

5  A(J,K)  =  A(J.K)/ROW(J) 

6  CONTINUE 

7  CONTINUE 

DETERM  =  CMPLXd.  0,0.0) 

SUMAXA  =  0. 0 
DO  20  J  =  1,  N 

ALPHA(J)  =0.0 
SUMROW  =  0.  0 
DO  10  I  =  1,  N 

ALPHA(J)  =  ALPHA(J)  +  A(J,I)*  CONJG(A( J,I)) 

10  SUMROW  =  SUMROW  +  CABS(A(J,I)) 

ALPHA'!)  =  SQRT(ALPHA(J)) 

IFCSUMROW.GT.  SUMAXA)  SUMAXA  =  SUMROW 
20  IPIVOT(J)  =  0 

DO  600  I  =  1,  N 
C 
C 

AMAX  =  CMPLX(0.  0,0.  0) 

DO  105  J  =  1,  N 

IF  (IPIVOT(J)-l)  60,  105,  60 
60  DO  100  K  =  1,  N 

IF  (IPIVOT(K)-l)  80,  100,  740 

80  TEMP  =  AMAX*CONJG(AMAX)  -  A( J,K)*CONJG(A( J,K)) 

IF( TEMP) 85, 85, 100 

85  IROW  =  J 

ICOLUM  s  K 
AMAX  =  A(J,K) 

100  CONTINUE 

105  CONTINUE 

IPIVOTC ICOLUM)  *  IPIVOTC ICOLUM)  +  1 
C 
C 

IF  (IROW-ICOLUM)  140,  260,  140 
140  DETERM  =  -DETERM 

DO  200  L  =  1,  N 

SWAP  =  ACIROW.L) 

ACIROW.L)  =  A(ICOLUM,L) 

200  A( ICOLUM, L)  =  SWAP 

SWAP  =  ALPHA(IROW) 

ALPHA(IROW)  =  ALPHAC ICOLUM) 

ALPHA( ICOLUM)  =  SWAP 
260  INDEX(I,1)  =  IROW 

INDEX(I,2)  =  ICOLUM 
PIVOT(I)  =  A( ICOLUM, ICOLUM) 

U  =  PIVOT(I) 

DETERM  =  DETERM*U 

DETERM  =  DETERM/ALPHA( ICOLUM) 

TEMP  =  PIVOT(I)*CONJG(PIVOT(I)) 

IF(TEMP)330, 720,330 
C 
C 

330  A( ICOLUM, ICOLUM)  =  CMPLX( 1. 0 ,0.  0) 
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350 

C 

DO  350  L  «  1,  N 

U  =  PIVOT(I) 

A(ICOLUM,L)  =  A(ICOLUM,L)/U 

* 

C 

380 

DO  550  LI  =  1,  N 

400 

IF(H-ICOLUM)  400,  550,  400 

T  =  A(L1,IC0LUM) 

A(Ll.ICOLUM)  =  CMPLX(0.0,0.  0) 

DO  450  L  =  1,  N 

450 

U  =  A(ICOLUM,L) 

A(L1,L)  =  A(L1,L)  -  U*T 

550 

CONTINUE 

600 

C 

CONTINUE 

C 

620 

DO  710  I  =  1,  N 

630 

L  =  N  +  1  -  I 

IF  (INDEX(L,1)  -  INDEX(L,2))  630,  710,  630 

JROW  =  INDEX(L,1) 

705 

JCOLUM  =  INDEX(L,2) 

DO  705  K  =  1,  N 

SWAP  =  A(K,JROW) 

ACK.JROW)  =  A(K,JCOLUM) 

A(K, JCOLUM)  =  SWAP 

CONTINUE 

710 

CONTINUE 

SUMAXI  =  0.0 

DO  910  I  =  1,  N 

SUMROW  =  0.  0 

900 

DO  900  J  =  1,  N 

SUMROW  =  SUMROW  +  CABS(A(I,J)) 

m 

910 

IFCSUMROW.GT.  SUMAXI)  SUMAXI  =  SUMROW 

CONTINUE 

950 

COND  =  SUMAXA*SUMAXI 

IFdNORM.NE.  1)  GO  TO  955 

DO  950  K  =  1,  N 

DO  950  J  =  1,  N 

A(J,K)  =  A(J,K)/(ROW(K)*COL(J)) 

955 

CONTINUE 

720 

RETURN 

WRITE(*,730) 

FORMATC  '  MATRIX  IS  SINGULAR') 

730 

740 

RETURN 

C 

C 

END 

SUBROUTINE  SAVE ( BCOND , ANS , U , OFFSET , PERND , CHAR5 , DPER , K , XORG , YORG 
CNRES , MRES , LBIAS , GBIAS , MAXD) 

THIS  SUBROUTINE  SAVES  THE  ESSENCE  OF  THE  FINITE  ELEMENT  PROBLEM 

C 

TO  A  DATA  FILE  CALLED  "  F3.DAT  tHIS  DATA  IS  NECESSARY  TO 

- 

C 

p 

SOLVE  THE  FIELD  FEEDBACK  FORMULATION,  (F3). 

COMPLEX  BCOND( 100) ,ANS( 100) ,U( 100,100) 

REAL  OFFSET , MESH( 0 :  2  00 , 5 ) , DPER , K , XORG , YORG , MRES , MAXD 

INTEGER  PERND,I,J,NRES,LBIAS,GBIAS 

119 


c 


CHARACTER*!  CHARS 


COMMON/BLKl/MESH 

C 

IF( (CHARS.  EQ.  'M' ).  OR. (CHARS. EQ.  ’m’))  THEN 
RETURN 

ENDIF 

C 

WRITE (40,*)  PERND 
WRITE (40,*)  OFFSET 
WRITE (40,*)  DPER 
WRITE(40,*)  K 
WRITE (40,*)  XORG 
WRITE(40,*)  YORG 
WRITE(40,*)  NRES 
WRITE(40,*)  MRES 
WRITE(40,*)  LBIAS 
WRITE(40,*)  GBIAS 
WRITE(40,*)  MAXD 
C 

DO  10,  I  =  1,  4 

DO  10,  J  =  1,  PERND 
WRITE(40,*)  MESH(J,I) 

10  CONTINUE 
C 

DO  20,  J  =  1,  PERND 

WRITE(40,*)  BCOND(J) 

20  CONTINUE 

C 

DO  30,  J  =  1,  PERND 

WRITE(40,*)  ANS(J) 

30  CONTINUE 

C 

DO  40,  1=1,  PERND 

DO  40,  J  =  1,  PERND 

WRITE(40,*)  U(I,J) 

40  CONTINUE 

C 

RETURN 

END 

C 

C 
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APPENDIX  D.  VARINT  CONVERGENCE  PROGRAM 


TEST  OF  VARINT  CONVERGENCE 

COMPLEX  F(3,3),  ALPHA,  BETA,  EXACT,  CENTER,  LEFT,  RIGHT,  TOP 

COMPLEX  BOTTOM,  SUM,  CALC 

REAL  X(3),  Y(3),  D,  AREA,  KR,  PI,  ERROR 

INTEGER  I 

OPEN  (UNIT  =  1,  FILE  =  'C:  MSFORT  TEST.DAT',  STATUS  =  'UNKNOWN’) 
ALPHA  =  (1. 0,0. 0) 

BETA  =  (1.0,0.  0) 

PI  =  4. 0*ATANC1.  0) 

DO  5,  I  =  1,  100 

D  =  2.  0*PI*FLOAT(I)/100.  0 
KR  =  REAL(CSQRT(BETA)) 

X(l)  =  0.0 
Y(l)  =0.0 
X(2)  =  D 
Y(2)  =  0.0 
X(3)  =  0.0 
Y(3)  =  D 

CALL  VARINT(X,Y,F, ALPHA, BETA, AREA) 

EXACT  =  CMPLX(C0S(KR*0.0),SIN(KR*0,0)) 

CENTER  =  4. 0*F(1,1) 

RIGHT  =  -2.  0*F(l,2)*CMPLX(C0S(KRn(2)),SIN(KR*Y(2))) 

TOP  =  -2.  0*F(1,2)*CMPLX(COS(KR*Y(3)),SIN(KR*Y(3))) 

LEFT  =  -2.  0*F(1,2)*CMPLX(COS(KR*Y(2)),SIN(KR*Y(2))) 

BOTTOM  =  -2.  0*F(1,2)*CMPLX(COS(-KR*Y(3)),SIN(-KR*Y(3))) 

SUM  =  TOP  +  BOTTOM  +  LEFT  +  RIGHT 
CALC  =  SUM/ CENTER 
ERROR  =  (CALC -EXACT) /EXACT 
WRITE( 1,1000)  I, D, EXACT, CALC, ERROR 
CONTINUE 

CLOSE(l) 

000  F0RMAT( IX, 13, IX, F8.  5 , 1X,2(F8.  5, 1X,F8.  5 , IX) ,E13.  6) 

STOP 

END 


SUBROUTINE  VARINT( X , Y , F , ALPHA , BETA , AREA ) 

GENERATING  VARIATIONAL  FINITE  ELEMENT  AREA  INTEGRATIONS  OF  THE 
LINEAR  BASIS  FUNCTION  LAGRANGIAN  FOR  THE  HELMHOLTZ  EQUATION. 
THESE  ARE  RETURNED  IN  F(3,3).  X(3)  AND  Y(3)  ARE  THE  WAVENUMBER 
NORMALIZED  CARTESIAN  COORDINATES  OF  THE  TRIANGLE  VERTICES. 

X  =  Ko*x,  Y  =  Ko*y  ALPHA  AND  BETA  ARE  COMPLEX 
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MATERIAL  PARAMETERS  WITHIN  THE  ELEMENT. 


FOR  TM  INCIDENCE:  ALPHA  =  1/ur;  BETA  =  er 
FOR  TE  INCIDENCE:  ALPHA  =  1/er;  BETA  =  ur 

OUTPUTS  :  F(3,3)  -  FINITE  ELEMENT  AREA  INTEGRATION 

AREA  -  AREA  OF  A  TRIANGLE 

COMPLEX  ALPHA,  BETA,  B12,  F(3,3) 

REAL  X(3),  y(3),  T(3,3),  AREA,  DET 
INTEGER  L,  K 

DET  =  ABSCX(2)*y(3)  +  X(3)*y(l)  +  X(l)*y(2)  -  X(3)*y(2)  - 
cx(i)*y(3)  -  x(2)*y(i)) 

AREA  =  ABS(0.  5*DET) 

B12  =  BETA/ 12. 

T(l,l)  =  (y(2)  -  y(3))/DET 

T(l,2)  =  (y(3)  -  y(l))/DET 

T(l,3)  =  (y(l)  -  y(2))/DET 

T(2,l)  =  (X(3)  -  X(2))/DET 

T(2,2)  =  (X(l)  -  X(3))/DET 

T(2,3)  =  (X(2)  -  X(1))/DET 

T(3,l)  =  (X(2)*y(3)  -  X(3)*y(2))/DET 

T(3,2)  =  (X(3)*y(l)  -  X(l)*y(3))/DET 

T(3,3)  =  (X(l)*y(2)  -  X(2)*y(l))/DET 

DO  10,  K  =  1,  3 

DO  10  L  =  1,  3 

F(K,L)  =  ALPHA*(T(1,K)*T(1,L)  +  T(2,K)*T(2,L))  -  B12 
IFCK.EQ.  L)  F(K,L)  =  F(K,L)  -  B12 
0  F(K,L)  =  AREA*F(K,L) 

RETURN 

END 
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APPENDIX  E.  FIELD  FEEDBACK  PROGRAM 


FIELD  FEEDBACK  FORMULATION  PROGRAM 
WRITTEN  BY  T.  B.  WELCH 

w/  PROGRAMMING  IDEAS  FROM  PROF  M. A.  MORGAN 
BCOND  -  BOUNDARY  CONDITIONS 

OFFSET  -  OFFSET  IN  WAVELENGTHS  (PERIMETER  TO  BOUNDARY) 

ANS  -  CALCULATED  PS I  VALUES  ON  PERIMETER 

DPER  -  DESIRED  PERCENT  ERROR  SCALE  FACTOR  FOR  GREEN'S 

FUNCTION  INTEGRAL  PATCH  STEPPING 
LBIAS  -  BIAS  THAT  IS  ADDED  TO  THE  GFI  STEP  FOR  <  1.  0 

GBIAS  -  BIAS  THAT  IS  ADDED  TO  THE  GFI  STEP  FOR  >1.0 

MAXD  -  MAXIMUM  DISTANCE  BEYOND  WHICH  NO  CONTRIBUTION  IS  MADE 

TO  THE  GFI 
K  -  WAVENUMBER 

XORG  -  X  ORIGIN 

YORG  -  Y  ORIGIN 

NRES  -  NUMBER  OF  EVENLY  SPACED  POINTS  DESIRED  FOR  THE  FAR 

FIELD  CALCULATIONS  (360/NRES  =  ANGULAR  RESOLUTION) 

U  -  MATRIX  THAT  RELATES  EACH  BOUNDARY  NODE  VALUE  TO  THE 

UNKNOWN  PERIMETER  NODE  VALUE.  MULTIPLY  U  BY  A  DRIVING 
VECTOR  (ON  BOUNDARY)  TO  FIND  PERIMETER  VALUES. 

PERND  -  NUMBER  OF  PERIMETER  NODES 

MESH  -  GEOMETRY  ARRAY  CONTAINING: 

1  -  X  POSITION  OF  PERIMETER  NODES 

2  -  Y  POSITION  OF  PERIMETER  NODES 

3  -  X  UNIT  NORMAL  OF  PERIMETER  NODES 

4  -  Y  UNIT  NORMAL  OF  PERIMETER  NODES 

5  -  X  POSITION  OF  GEOMETRIC  CONTOUR  NODES 

6  -  Y  POSITION  OF  GEOMETRIC  CONTOUR  NODES 

7  -  X  POSITION  OF  BOUNDARY  NODES 

8  -  Y  POSITION  OF  BOUNDARY  NODES 

T  -  MATRIX  THAT  RELATES  PERIMETER  VALUES  BACK  OUT 

TO  THE  BOUNDARY  VIA  A  GREEN'S  FUNCTION  INTEGRAL 
CNVEC  -  VECTOR  OF  SCATTERED  FIELD  BACK  ONTO  THE  BOUNDARY 
MRES  -  MESH  RESOLUTION 


COMPLEX  BCOND( 100) ,ANS( 100) ,U( 100,100) ,T( 100 , 100) ,CNVEC( 100) 

REAL  OFFSET, MESH( 100,8) , DPER, K, XORG, YORG, MRES, MAXD 
INTEGER  PERND, I, J, NRES, LB IAS, GBIAS 
C 

OPEN  (UNIT  =  40,  FILE  =  'C;  MSFORT  F 3. DAT' ,STATUS=' UNKNOWN ' ) 

OPEN  (UNIT  =  50,  FILE  =  'C:  MSFORT  FFPAT.  DAT ', STATUS= ' UNKNOWN ' ) 

C 

CALL  INPUT( BCOND , ANS , U , OFFSET , PERND , MESH , DPER , K , NRES , XORG , YORG , 
CMRES , LBIAS , GBIAS , MAXD) 

CALL  TMAT( U , PERND , MESH , T , OFFSET , DPER , BCOND , MRE S , LB IAS , GB IAS , MAXD ) 
CALL  CNSOLV(T, BCOND, CNVEC, PERND) 
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CALL  FFLD( CNVEC , PERND , MESH , U , OFFSET , K , NRES , XORG , YORG ) 
C 

CLOSE (40) 

CL0SE(50) 

C 

STOP 

END 


SUBROUTINE  INPUT( BCOND , ANS ,U , OFFSET , PERND , MESH , DPER , K , NRES , XORG , 
CYORG , MRES , LB IAS , GB IAS . MAXD ) 

THIS  SUBROUTINE  READS  THE  FINITE  ELEMENT  PROBLEM  DATA  FROM 
THE  DATA  FILE  CALLED  "  F3.DAT  THIS  DATA  IS  NECESSARY  TO 
SOLVE  THE  FIELD  FEEDBACK  FORMULATION,  (F3). 


COMPLEX  BCOND( 100) , ANS( 100) ,U( 100,100) 

REAL  OFFSET , MESH( 1 00 , 8 ) , DPER ,K ,XORG , YORG , MRES , MAXD 
INTEGER  PERND,I,J,NRES,LBIAS,GBIAS 


WRITEC*,*) 

'  READING  INPUT  DATA  ' 

1 

READ(40,*) 

PERND 

READ(40,*) 

OFFSET 

READ(40,*) 

DPER 

READ(40,*) 

K 

READ(40,*) 

XORG 

READ(40,*) 

YORG 

READ(40,*) 

NRES 

READ(40,*) 

MRES 

READ(40,*) 

LB  IAS 

READ(40,*) 

GBIAS 

READ(40,*) 

MAXD 

- 

WRITE(6,*) 

'  NUMBER  OF  PERIMETER  NODES 

= 

' , PERND 

WRITE(6,*) 

'  BOUNDARY  CONTOUR  OFFSET 

S 

' .OFFSET 

WRITE(6,*) 

'  DESIRED  GFI  SCALE  FACTOR 

= 

' ,DPER 

WRITE(6,*) 

'  GFI  STEP  BIAS  FOR  <1.0 

s 

’ , LB IAS 

WRITE(6,*) 

'  GFI  STEP  BIAS  FOR  >  1.  0 

s 

' .GBIAS 

WRITE(6,*) 

'  MAX  DIST  >  NO  CONTRIBUTION  TO  GFI 

= 

' .MAXD 

VRITE(6,*) 

'  WAVENUMBER 

s 

',K 

WRITEC 6,*) 

'  X  ORIGIN 

= 

' ,XORG 

WRITE(6,*) 

'  Y  ORIGIN 

s 

' .YORG 

WRITEC 6,*) 

'  NUMBER  OF  NODES  FOR  SIGMA 

= 

' .NRES 

WRITEC6,*) 

'  REQUESTED  MESH  RESOLUTION 

s 

' .MRES 

DO  10,  I  = 

1,  A 

DO  10 

,  J  =  1,  PERND 

READ(40,*)  MESH(J,I) 

0  CONTINUE 

DO  20,  J  =  1,  PERND 

READ(40,*)  BCOND(J) 

20  CONTINUE 

C 

DO  30,  J  =  1,  PERND 
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READ(40,*)  ANS(J) 

30  CONTINUE 

C 

DO  40,  1=1,  PERND 

DO  40,  J  =  1,  PERND 

READ(40,*)  U(I,J) 

40  CONTINUE 

C 

DO  50,  1=1,  PERND 

MESH(I,5)  =  MESH(I,1)  +  MESH( I ,3)*OFFSET/2. 0 
MESH(I,6)  =  MESH(I,2)  +  MESH( 1 ,4)*0FFSET/2.  0 
MESH(I,7)  =  MESH(I,1)  +  MESH( I ,3)*0FFSET 
MESH(I,8)  =  MESH(I,2)  +  MESH( I ,4)*0FFSET 
50  CONTINUE 

C 

RETURN 

END 

C 

C 

SUBROUTINE  TMAT( U , PERND , MESH , T , OFFSET , DPER , BCOND , MRES , LB I AS , GB I AS , 
CMAXD) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  GREEN’S  FUNCTION  INTEGRAL 
C  (FOR  A  SINGLE  BASIS  FUNCTION  BOUNDARY  CONDITION)  GEOMETRIC 

C  PERIMETER  WITH  RESPECT  TO  EACH  OF  THE  OFFSET  BOUNDARY  NODES. 

C  THE  INTEGRATION  IS  REPEATED  UNTIL  EACH  BOUNDARY  CONDITION  HAS  BEEN 
C  INDIVIDUALLY  APPLIED  AND  INTEGRATED  WITH  RESPECT  TO  EACH  OFFSET 
C  BOUNDARY  NODE.  THESE  VALUES  ARE  RETURNED  IN  THE  ”  T  "  MATRIX 
C  FOR  USE  IN  THE  FIELD  FEEDBACK  FORMULATION.  THE  MATRIX  IS 
C  ORGANIZED,  T  m,n. 

C 

COMPLEX  U( 100,100) ,T( 100,100) ,PVEC( 100) ,PDVEC( 100) 

COMPLEX  J,H0RP1,H0RP2,H1RP1,SUM,TEMP(100),BC0ND(100) 

COMPLEX  H1RP2 , PS I , PS IRP , INTEGRAL , DPS I C 
REAL  RMRP,N0RM11,N0RM12,D0T,DPER,DISTM,MAXD 
REAL  OFFSET, MESH( 100,8) ,DIST,R,DZ,DL,MRES 

INTEGER  I,M,N,NN,STEP,FN,SN,PERND,MM,STEPMX,STEPMN,LBIAS,GBIAS 
OPEN  (UNIT  =  2,  FILE  =  ’C:  MSFORT  TMAT.  DAT' , STATUS  =  'UNKNOWN') 

C 

J  =  (0.0, 1.0) 

STEPMX  =  INT(DPER*( -48.  0*0FFSET  +  17.2)  +  LBIAS) 

STEPMN  =  INT(DPER  +  GBIAS) 

C 

WRITE(*,*)  '  LOADING  T  MATRIX  ' 

WRITE(*,*)  '  MAXIMUM  STEP  =  ', STEPMX,’,  MINIMUM  STEP  =  ', STEPMN 
WRITE(*,*)  '  MAXIMUM  DISTANCE  FOR  ANY  CONTRIBUTION  =  ’ ,MAXD 
DO  40,  M  =  1,  PERND 

DO  5,  1=1,  PERND 
IF(M.EQ.  I)  THEN 

PVEC(I)  =  (1.0  +  U(I,M))/2.0 
PDVEC(I)  =  (1.0  -  U( I, M)) /OFFSET 

ELSE 

PVEC(I)  =  U(I,M)/2. 0 
PDVEC(I)  =  -U( I, M) /OFFSET 
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ENDIF 

5  CONTINUE 

C 

DO  30,  N  »  1,  PERND 

WRITE(*,1000)  M,  N,  PERND 
SUM  «  (0.  0,0.  0) 

DO  20,  NN  »  1,  PERND 
FN  »  NN 

IFCNN.EQ. PERND)  THEN 
SN  »  1 

ELSE 

SN  «  NN  +  1 

ENDIF 

DIST  »  SQRT((MESH(FN,5)-MESH(SN,5))**2+(MESH(FN,6) 

C  -MESH(SN,6))**2) 

INTEGRAL  »  (0.  0,0.  0) 

C 

DISTM  =  SQRT((MESH(N,7)-MESH(FN,5))**2+(MESH(N,8)-MESH(FN,A))*^2) 

C 

R  =  MESH(SN,5)  -  MESH(FN,5) 

DZ  =  MESH(SN,6)  -  MESH(FN,6) 

DL  =  SQRT(R**2  +  DZ**2) 

NORMll  =  -DZ/DL 
N0RM12  =  R/DL 
C 

IF(DISTM.GT.MAXD)  THEN 
GOTO  20 

ELSEIF(DISTM.  LE.  1.  0)  THEN 
STEP  =  STEPMX 

ELSE 

STEP  =  STEPMN 

ENDIF 

IF(STEP.LT.  1)  STEP  =  1 
C 

DO  10,  I  =  1,  STEP+1 
IFCI.EQ.  1)  THEN 

RMRP  =  SQRT((MESH(N,7)-(MESH(FN,5)+0.  25*(MESH(SN,5)- 
MESH( FN , 5 ) ) /FLOAT( STEP) ) )**2+( MESH( N , 8 ) - ( MESH( FN , 6 ) + 

0.  25*( MESH( SN , 6 ) -MESH( FN , 6 ) ) /FLOAT( STEP) ) )**2 ) 

DPSIC  =  PDVEC(FN)  +  0.  25*(PDVEC(SN)-PDVEC(FN))/STEP 
PSIRP  =  PVEC(FN)  +  0. 25*(PVECCSN)-PVEC(FN))/STEP 
DOT  =  (N0RM11*(MESH(N,7)  -  (MESH(FN,5)+0.25*(MESH(SN,5)- 
MESH(FN,5))/STEP))+(NORM12*(MESH(N,8)-(MESH(FN,6)+0.  25* 

( MESH( SN , 6 ) -MESH( FN , 6 ) ) /STEP ) ) ) ) /RMRP 
ELSEIFC I. EQ.  STEP+1)  THEN 

RMRP  =  SqRT((MESH(N,7)-(MESH(FN,5)+(FL0AT(STEP)-0.  25)*( 
MESH( SN , 5 ) -MESH( FN , 5 ) ) /FLOAT( STEP) ) ) **2+( MESH( N , 8 ) - ( MESH 
( FN , 6 )+( FLOAT( STEP) -0.  25 )*( MESH( SN , 6 ) -MESH( FN , 6 ) ) /FLOAT 
(STEP)))**2) 

DPSI C=PDVEC ( FN ) +( FLOAT( STEP ) - 0 . 25 ) *( PDVEC( SN ) - PDVEC ( FN ) ) 
/(STEP) 

PSIRP  =  PVEC(FN)+(FLOAT(STEP)-0. 25)*(PVEC(SN)-PVEC(FN))/ 
(STEP) 

DOT  =  (NORM11*(MESH(N,7)-(MESH(FN,5)+(FLOAT(STEP)-0.  25)* 
( MESH( SN , 5 ) -MESH(FN , 5 ) ) /STEP) )+( N0RM12*( MESH( N , 8 ) - ( MESH 
(FN,6)+(FL0AT(STEP)-0.  25)*(MESH(SN,6)-MESH(FN,6))/STEP) 
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ELSE 

RMRP  =  SQRT((MESH(N,7)-(MESH(FN.5)+FLOAT(I-l)*(MESH(SN. 

C  5 ) *MESH( FN , 5 ) ) /FLOAT( STEP) ) )**2+( MESH( N , 8 ) -( MESH( FN . 6 )+ 

C  FL0AT(I-1)*(MESH(SN,6)-MESH(FN.6))/FL0AT(STEP)))**2) 

DPS IC*PDVEC( FN )+( I - 1 )*( PDVEC( SN ) -PDVXC( FN ) ) / ( STEP ) 

PSIRP  =  PVEC(FN)+(I-1)*(PVEC(SN)-PVEC(FN))/(STEP) 

DOT  =  (N0RM11*(MESH(N,7)-(MESH(FN,5)+(I-1)*(MESH(SN,5)- 
C  MESH(FN,5))/STEP))+(N0RM12*(MESH(N,8)*(KESH(FN,6)+(I-1)* 

C  ( ffiSH( SN , 6 ) -MESH( FN . 6 ) ) /STEP) ) ) ) /RMRP 

ENDIF 

CALL  HAN1(CMPLX(RMRP,0.  0),HORP1,H1RP1) 

IF((I.EQ.  l).OR.  (I.EQ.  STEP+D)  THEN 

PSI  =  (J/4. 0)*(H0RP1*DPSIC  -  PSIRP*DOT*HlRPl)/2 

ELSE 

PSI  »  (J/4.  0)*(H0RP1*DPSIC  -  PSIRP*D0T*H1RP1) 

ENDIF 

INTEGRAL  =  INTEGRAL  +  PSI*DIST/STEP 
C 

10  CONTINUE 

SUM  =  SUM  +  INTEGRAL 
20  CONTINUE 

T(N,M)  =  SUM 
30  CONTINUE 

40  CONTINUE 
C 

DO  55,  1=1,  PERND 

TEMP(I)  =  (0. 0,0.0) 

DO  50,  M  =  1,  PERND 

TEMP(I)  =  TEMP(I)  +  T(I,M)*BCOND(M) 

50  CONTINUE 

WRITE(2,*)  (T(I,MM),  MM  =  1,  PERND) 

55  CONTINUE 

C 

DO  60,  I  =  1,  PERND 

WRITE(2,*)  I,  TEMP(I) 

60  CONTINUE 

C 

1000  FORMATC IX, 'COLUMN  ',13,',  ROW  ',13,'  OUT  OF  ',13) 

C 

CL0SE(2) 

C 

RETURN 

END 

C 

C 

SUBROUTINE  CNSOLV(T,BCOND,CNVEC, PERND) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  C  VECTOR  BY  SOLVING: 

C 

C  -1 

C  Cn  =  [I  -  T]  *  BOUNDARY  CONDITIONS  (INCIDENT  FIELDS) 

C 

COMPLEX  BCOND(100),T(100,100),TEMP(100),CNVEC(100),  DETERM 
REAL  COND,  DMAG 


n  o  o  o  o  Cl  o 


INTEGER  PERND,I,J,K,L,INORM,NMAX 
C 

INORM  =  0 
NMAX  =  100 
C 

DO  10,  I  =  1,  PERND 

DO  10,  J  =  1,  PERND 
IF(I.EQ.J)  THEN 

T(I,J)  =  1  -  T(I,J) 

ELSE 

T(I,J)  =  -T(I,J) 

ENDIF 

10  CONTINUE 

WRITE(*,*)  'INVERTING  THE  [  I  -  T]  MATRIX 
C 

CALL  CSMINVC T , NMAX , PERND , DETERM , COND , INORM) 

C 

DMAG  =  CABS (DETERM) 

WRITEC*,*)  '  DETERMINANT  =  ',  DMAG 

WRITEC*,*)  '  CONDITION  NUMBER  =  ’,  COND 

WRITEC*,*)  '  MULTIPLING  MATRICES  TO  FORM  THE  SCATTERED  FIELDS  ' 
C 

DO  30,  1=1,  PERND 

CNVEC(I)  =  (0. 0,0.0) 

DO  30,  J  =  1,  PERND 

CNVEC(I)  =  CNVEC(I)  +  T( I , J)*BCOND( J) 

30  CONTINUE 
C 

RETURN 

END 


SUBROUTINE  FFLD( CNVEC , PERND , MESH , U , OFFSET , K , NRES , XORG , YORG ) 

THIS  SUBROUTINE  CALCULATES  THE  FAR  FIELDS  DUE  TO  THE  OFFSET 
BOUNDARY  SCATTERED  FIELDS  AND  THE  PERIMETER  SCATTERED  FIELDS. 
ADDITIONAL  GREEN'S  FUNCTION  INTEGRALS  ARE  ACCOMPLISHED. 

COMPLEX  CNVEC( 100) ,U( 100,100) , J,PSISP( 100) ,TEMP,PSI ,DPSI 
COMPLEX  DPSISPC 100) , INTEGRAL, DPSIM( 100 , 100) ,PSIM( 100 , 100) 

REAL  MESHC 100 , 8 ) , OFFSET , K , DOT , DOTl , PI , ARES , XORG , YORG , D I ST , S IGMA 
REAL  R,DZ,DL,N0RM11,N0RM12 
INTEGER  PERND, I, M,N,L, NRES, FN,SN, STEP, II 
C 

WRITEC*,*)  '  CALCULATING  THE  SCATTERED  FIELDS' 

J  =  (0.  0,1.  0) 

PI  =  4.  0*ATAN(1. 0) 

ARES  =  2.  0*PI/FL0AT(NRES) 

STEP  =  5 
C 

DO  5,  M  =  1,  PERND 

DO  5,  I  =  1,  PERND 
IFCM.EQ. I)  THEN 

PSIM(I,M)  =  (1.0  +  U(I,M))/2.0 
DPSIM(I,M)  =  (1.0  -  U(I,M))/OFFSET 

ELSE 
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PSIMd.M)  =  U(I,M)/2.  0 
DPSIM(I.M)  =  -U( I, M) /OFFSET 

ENDIF 

5  CONTINUE 
C 

DO  10.  1=1,  PERND 

PSISP(I)  =  (0. 0,0. 0) 

DPSISP(I)  =  (0. 0,0. 0) 

DO  10,  L  =  1,  PERND 

PSISP(I)  =  PSISP(I)  +  PSIM(I,L)*CNVEC(L) 

DPSISP(I)  =  DPSISP(I)  +  DPSIM(I,L)*CNVEC(L) 

10  CONTINUE 

C 

DO  30,  1=0,  NRES-1 

WRITE(6,1000)  I+l,  NRES 
INTEGRAL  =  (0. 0,0. 0) 

DO  20,  N  =  1,  PERND 
FN  =  N 

IFCN.EQ.  PERND)  THEN 
SN  =  1 

ELSE 

SN  =  N  +  1 

ENDIF 

R  =  MESH(SN,5)  -  MESH(FN,5) 

DZ  =  MESH(SN,6)  -  MESH(FN,6) 

DL  =  SQRT(R**2  +  DZ**2) 

NORMll  =  -DZ/DL 

NORM12  =  R/DL 

DO  15,  II  =  1,  STEP+1 

DIST  =  SQRT((MESH(FN,5)-MESH(SN,5))**2+(MESH(FN,6)- 
C  MESH(SN,6))’«rt^2) 

IFdI.EQ.  1)  THEN 

PSI  =  PSISP(FN)+0. 25*(PSISP(SN)-PSISP(FN))/ 

C  FLOAT(STEP) 

DPSI  =  DPSISP(FN)+0.  25*(DPSISP(SN)*DPSISP 
C  (FN))/FLOAT(STEP) 

DOT  =  NORMll*SINd*ARES)  +  NORM12*COSd*ARES) 

DOTl  =  (MESH(FN,5)+0. 25*(MESH(SN,5)-MESH(FN, 
5))/FLOAT(STEP)-XORG)*SINd*ARES)  +  (MESH(FN,6)-t- 
0.  25*( MESH( SN, 6 ) -MESH( FN , 6) ) /FLOAT( STEP) - 
YORG)*COSd*ARES) 

ELSEIFCII.EQ. STEP+1)  THEN 

PSI  =  PSISP(FN)+(FL0AT(STEP)-0.25)*(PSISP(SN)- 
PS I SP( FN ) ) /FLOAT( STEP ) 

DPSI  =  DPSISP(FN)+(FL0AT(STEP)-0.25)*(DPSISP(SN)- 
DPSISP(FN) ) /FLOAT( STEP) 

DOT  =  NORMll*SINd*ARES)  +  N0RM12*C0Sd*ARES) 

DOTl  =  (MESH(FN,5)+(FLOAT(STEP)-0.  25)*(MESH(SN,5)- 
MESH(FN,5))/FLOAT(STEP)-XORG)*SINd*ARES)  +  (MESH( 
FN,6)+(FL0AT(STEP)-0.  25)*(MESH(SN,6)-MESH(FN,6))/ 
FLOAT( STEP) -YORG)*COS( I*ARES) 

ELSE 

PSI  =  PSISP(FN)+FLOATdI-l)*(PSISP(SN)-PSISP(FN))/ 
C  FLOAT(STEP) 

DPSI  =  DPSISP(FN)+FLOATdI-l)*(DPSISP(SN)-DPSISP 
C  (FN))/FLOAT(STEP) 
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DOT  =  NORMll*SIN(I*ARES)  +  NORM12*COS( I*ARES) 

DOTl  =  (MESH(FN,5)+FLOAT(II-l)*(MESH(SN,5)-MESH(FN. 
C  5))/FLOAT(STEP)-XORG)*SIN(I*ARES)  +  (MESH(FN,6)+ 

C  FLOAT(II-l)*(MESH(SN.6)-MESH(FN,6))/FLOAT(STEP)- 

C  YORG)*COS(I*ARES) 

ENDIF 

IF((II.EQ.  l).OR.(II.EQ.  STEP+D)  THEN 

TEMP=(J*DPSI+DOT*PSI)*(EXP(J*DOTl))*DIST/(2.  0* 
C  STEP) 

ELSE 

TEMP=( J*DPSI+DOT*PSI )*( EXP( J*D0T1 ) )*DIST/ 

C  (STEP) 

ENDIF 

INTEGRAL  =  INTEGRAL  +  TEMP 
15  CONTINUE 

20  CONTINUE 

SIGMA  =  ((CABS(INTEGRAL))**2.  0)/(4.  0*K) 

WRITE(50,1010)  I+l,  (I*ARES*180. 0/PI),  SIGMA 
30  CONTINUE 
C 

WRITE(6,*)  ’ 

WRITE(6,*)  '  <«  FIELD  PATTERN  STORED  IN  FFPAT.DAT  »>  ’ 

VnRITE(6.*)  '  ' 

C 

1000  FORMATC IX, 'INTEGRAL  '.13,',  OUT  OF  ’,13,'  COMPLETED') 

1010  F0RMAT(1X,I3,2X,F6.  2,2X,E14.  8) 

C 

RETURN 

END 

C 

C 

SUBROUTINE  HANK Z, HO, HI) 

C 

C  Computing  Hankel  Functions  for  n=0,l  with 
C  Complex  Argument,  Z.  Direct  Power  Series  Method  for 
C  CABS(Z)  .  LE.  5  and  Hankel's  Asymptotic  Formula  for 
C  CABS(Z)  .GT.  5.  Written  11/6/87  by  M. A.  Morgan 
C 

INTEGER  M,M2 

REAL  C(34) ,DM,F(34) ,G0,P(34) ,Pi,P2 

COMPLEX  Z,Z2,Z3,Z4,J0,J1,Y0,Y1,AM,CL,P0,P1,Q0,Q1 

COMPLEX  E0,El,X0,Xl,H0,Hl,j 

PI=3.  1415927 

P2=2. 0/PI 

j=(0.  .1.  ) 

IF(CABS(Z).LE.  5.0)  THEN 
C 

C  Direct  Power  Series  Method 

C 

G0=  1.78072 
Z2=0.  5*Z 
CL=CL0G(G0*Z2) 

C 

C  Computing  F(m)  =  m  !  and  P(m)  =  1  +  1/2  +  1/3  +  _ +  1/m 

C 

F(l)=l. 0 
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P(l)=l.  0 
DO  11  M=2,34 

F(M)=M*F(M-1) 

P(M)=P(M-1)+1.0/M 
11  CONTINUE 

C 

C  Computing  Power  Series  Coefficients 
C 

DM=-1.  0 
DO  22  M=l,34 

C(M)=DM/(F(M)*F(M)) 

DM=-DM 

22  CONTINUE 

C 

C  Computing  JO  and  J1 
C 

J0=(1.  ,0.  ) 

J1=(0.  ,0.  ) 

M=0 

33  M=M+1 

M2=2*M 

AM=C(M)*(Z2**M2) 

J0=J0+AM 

J1=J1-M*AM 

IF((CABS(AM).GT,  1.  0E-10).AND.  (M.  LT.  34))  GO  TO  33 
J1=J1/Z2 
C 

C  Computing  YO  and  Y1 
C 

M=0 

Y0=CL*J0 

Yl=Z2*CL*Jl-0.  5*J0 
44  M=M+1 

M2=2*M 

AM=C(M)*P(M)*(Z2**M2) 

Y0=Y0-AM 

Y1=Y1+M*AM 

IF((CABS(AM).GT.1.0E-10).AND.(M.LT.34))  GO  TO  44 

Y0=P2*Y0 

Y1=P2*Y1/Z2 

H0=J0-j*Y0 

Hl=Jl-j*Yl 

RETURN 

ELSE 

C 

C  Hankel'  Asymptotic  Formula  (Abram.  &  Stegun  p.  364 

c 

Z2=Z*Z 

Z3=Z*Z2 

Z4=Z*Z3 

P0=1. 0-.  0703125/Z2+.  1121521/Z4 
Q0=-.  125/Z+.  0732422/Z3 
Pl=l.  0+.  1171875/Z2-.  1441956/Z4 
Ql=.  375/Z-.  10253906/Z3 
X0=(Z-.  25*PI) 

X1=(Z-. 75*PI) 
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E0=CEXP(-j*X0) 

El=CEXP(-j*Xl) 

AM=CSQRT(P2/Z) 

H0=AM*(P0-j*Q0)*E0 

Hl=AM*(Pl-j*Ql)*El 

ENDIF 

C 

RETURN 

END 


SUBROUTINE  CSMINV(A,NDIM,N,DETERM,COND,INORM) 


INORM  -  FLAG  TO  NORMALIZE  COLUMNS  AND  ROWS  OF  MATRIX  A 
MATRIX  NORMALIZATION  BY  M.  A.  MORGAN 
APRIL  24,1978 


A  -  MATRIX  TO  INVERT 

NDIM 

N 

DETERM  -  DETERMINATE  OF  A 
COND  -  CONDITION  NUMBER  OF  A 
INORM  -  INTEGER  NORMALIZATION  FLAG 


INPUT/OUTPUT 

INPUT 

INPUT 

OUTPUT 

OUTPUT 

INPUT 


COMPLEX  A( 100, 100), PIVOTC 100), AMAX.T, SWAP, DETERM, U 
INTEGER  I ,  J , K ,  L ,  IPIV0T(  100 )  ,  INDEX(  100 , 2.)  ,  IROW ,  ICOLUM ,  LI ,  JROW 
INTEGER  JC0LUM,N, INORM 

REAL  TEMP , ALPHA( 100 ) , C0L( 100 ) , ROW( 100 ) , A JK , SUMAXA , SUMROW , SUMAXI 
IF(NDIM.  GT.  100)  THEN 

WRITE(*,*)  '  ERROR  IN  INVERTION  CALL...  DIMENSION  >  100  ’ 
STOP 

ENDIF 

IF(N. GT.NDIM)  THEN 

WRITEC*,*)  '  ERROR  IN  INVERTION  CALL...  N  >  MAX  DIM.  ' 

STOP 

ENDIF 

IF( INORM.  NE.  1)  GO  TO  7 
DO  3  K  =  1,  N 

COL(K)  =0.0 
DO  1  J  =  1,N 

AJK  =  CABS(A(J,K)) 

IF(AJK.GT.COL(K))  COL(K)  =  AJK 
CONTINUE 
DO  2  J  =  1,  N 
A(J,K)  =  A(J,K)/COL(K) 

CONTINUE 
ROW  NORMALIZING 
DO  6  J  =  1,  N 

ROW(J)  =0.0 
DO  4  K  =  1,  N 

AJK  =  CABS(A(J,K)) 

IF(AJK.GT.ROW(J))  ROW(J)  =  AJK 
4  CONTINUE 
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DO  5  K  *  1,  N 

5  A(J,K)  =  A(J,K)/ROW(J) 

6  CONTINUE 

7  CONTINUE 

DETERM  =  CMPLXCl.  0,0.0) 

SUMAXA  =  0.  0 
DO  20  J  =  1,  N 

ALPHA(J)  =  0.0 
SUMROW  =  0.  0 
DO  10  I  =  1,  N 

ALPHA(J)  =  ALPHA(J)  +  A(J,I)*  CONJG(A( J, I) ) 

10  SUMROW  =  SUMROW  +  CABS(A(J,I)) 

ALPHA(J)  =  5QRT(ALPHA(J)) 

IFCSUMROW.GT.  SUMAXA)  SUMAXA  =  SUMROW 
20  IPIVOT(J)  =  0 

DO  600  I  =  1,  N 
C 
C 

AMAX  =  CMPLX(0.  0,0.  0) 

DO  105  J  =  1,  N 

IF  (IPIVOT(J)-l)  60,  105,  60 
60  DO  100  K  =  1,  N 

IF  (IPIVOT(K)-l)  80,  100,  740 

80  TEMP  =  AMAX*CONJG(AMAX)  -  A( J,K)*CONJG(A( J,K)) 

IF( TEMP) 85, 85, 100 

85  IROW  =  J 

ICOLUM  =  K 
AMAX  =  A(J,K) 

100  CONTINUE 

105  CONTINUE 

IPIVOT(ICOLUM)  =  IPIVOT(ICOLUM)  +  1 

C 

IF  (IROW-ICOLUM)  140,  260,  140 
lAO  DETERM  =  -DETERM 

DO  200  L  =  1,  N 

SWAP  =  A(IROW,L) 

A(IROW,L)  =  A(ICOLUM,L) 

200  A(ICOLUM,L)  =  SWAP 

SWAP  =  ALPHA(IROW) 

ALPHA(IROW)  =  ALPHA(ICOLUM) 

ALPHA(ICOLUM)  =  SWAP 
260  INDEX(I,1)  =  IROW 

INDEX(I,2)  =  ICOLUM 
PIVOT(I)  =  A( ICOLUM, ICOLUM) 

U  =  PIVOT(I) 

DETERM  =  DETERM*U 

DETERM  =  DETERM/ ALPHA( ICOLUM) 

TEMP  =  PIVOT(I)*CONJG(PIVOTCI)) 

IF(TEMP)330,720,330 

C 

C 

330  A( ICOLUM, ICOLUM)  =  CMPLX( 1.0, 0.0) 

DO  350  L  =  1,  N 
U  =  PIVOT(I) 

350  A( ICOLUM, L)  =  A( ICOLUM, L)/U 
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c 

c 

380  DO  550  LI  =  1,  N 

IF(Ll-ICOLUM)  400,  550,  400 
400  T  =  A(L1,IC0LUM) 

A(L1,IC0LUM)  =  CMPLXCO.  0,0.  0) 

DO  450  L  =  1,  N 

U  =  A(ICOLUM,L) 

450  A(L1,L)  =  A(L1,L)  -  U*T 

550  CONTINUE 

600  CONTINUE 
C 
C 

620  DO  710  I  =  1,  N 

L  =  N  +  1  -  I 

IF  (INDEX(L,1)  -  INDEX(L,2))  630,  710,  630 
630  JROW  =  INDEX(L,1) 

JCOLUM  =  INDEX(L,2) 

DO  705  K  =  1,  N 

SWAP  =  A(K,JROW) 

A(K,JROW)  =  A(K, JCOLUM) 

A(K, JCOLUM)  =  SWAP 
705  CONTINUE 

710  CONTINUE 

SUMAXI  =  0.0 
DO  910  I  =  1,  N 
SUMROW  =  0.  0 

DO  900  J  =  1,  N 

900  SUMROW  =  SUMROW  +  CABS(A(I,J)) 

IFCSUMROW.GT.  SUMAXI)  SUMAXI  =  SUMROW 
910  CONTINUE 

COND  =  SUMAXA*SUMAXI 
IFCINORM.NE.  1)  GO  TO  955 
DO  950  K  =  1,  N 

DO  950  J  =  1,  N 

950  A(J,K)  =  A(J,K)/(ROW(K)*COL(J)) 

955  CONTINUE 
RETURN 

720  WRITE(*,730) 

730  FORMAT(  '  MATRIX  IS  SINGULAR') 

740  RETURN 
END 
C 
C 
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APPENDIX  F.  DIELECTRIC  CYLINDER  SCATTERING  PROGRAM 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

PLANEWAVE  SCATTERING  BY  A  DIELECTRIC  CYLINDER 
E  -  WAVE  (TM  CASE) 

H  -  WAVE  (TE  CASE) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

COMPLEX  GAMMA(0:  200,2) ,SIGMA(2) ,ER, MU, KR,YR,ZR 
COMPLEX*16  JA(0:200),DJA(0:200),KRRA 
C0MPLEX*16  J(0;  200) ,DJ(0: 200) 

REAL*8  Y(0: 200) ,DY(0:  200) ,YA(0:  200) ,DYA(0:  200) , JB(0;  200) 

REAL*8  DJB(0: 200), RA, KORA 
REAL  K0,PI,SIGMAN(2),A,B 
INTEGER  I, II, MODE, ARES, N 
C 

0PEN(3,FILE='C:  MSFORT  DECTEM.DAT') 

PI  =  3. 1415927 
C 

WRITE(*,*)  'PERMITTIVITY  FOPvMAT  IS  "  a  +  jb,  "  ' 

WRITE(*,*)  'Enter  Dielectric  Constant  (REAL  PART,  a)  ' 

READ(*,*)  A 

WRITEC*,*)  'Enter  Dielectric  Constant  (IMAGINARY  PART,  b)  ' 
READ(*,*)  B 
ER  =  CMPLX(A,B) 

C 

WRITE(*,*)  'PERMEABILITY  FORMAT  IS  "  a  +  jb,  "  ' 

WRITE(*,*)  'Enter  Permeability  Constant  (REAL  PART,  a)  ' 

READ(*,*)  A 

WRITE(*,*)  'Enter  Permeability  Constant  (IMAGINARY  PART,  b)  ' 
READ(*,*)  B 
MU  =  CMPLX(A,B) 

C 

WRITE(*,*)  'Enter  the  wave  number  (ko)  ' 

READ(*,*)  KO 
C 

WRITE(*,*)  'Enter  Cylinder  Radius  (IN  WAVENUMBER  UNITS)  ' 
WRITE(*,*)  'WARNING:  Do  Not  Enter  Zero  !  ' 

READ(*,*)  KORA 
C 

KR  =  CSQRT(MU*ER) 

YR  =  CSQRT(ER/MU) 

ZR  =  CSQRT(MU/ER) 

RA  =  KORA/KO 
KRRA  =  KR*RA 
C 

WRITE(*,*)  'Enter  No.  of  Modes;  ' 

READ(*,*)  MODE 

WRITE(*,*)  'Enter  the  angular  resolution  * 

READ(*,*)  ARES 
C 
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WRITE(3,*)  'Cylinder  Scattering  vs.  angle' 

WRITE(3,1101  ER,MU,RA,MODE,KO,KRRA 
WRITE(3,*)  ’  ' 

C 

CALL  BES(MODE,KORA,JB,Y,DJB,DY) 

CALL  DCBJNS  (DCMPLX(K0RA,0.  OD+00) .MODE, J,DJ) 

CALL  DCBJNS  ((KRRA*K0) .MODE.JA.DJA) 

WRITEC*,*)  '  RETURNED  FROM  FINAL  BESSEL  CALL' 

C 

C  CALCULATING  GAMMAs 

C 

DO  10,  N  =  0,  MODE 

GAMMA(N,1)  =  (JA(N)*DJ(N)-YR*DJA(N)*J(N))/(JA(N)*CMPLX(DJ(N) 
C  ,-DY(N))  -  YR*DJA(N)*CMPLX(J(N),-Y(N))) 

GAMMA(N,2)  =  ( JA(N)*DJ(N)-ZR*DJA(N)*J(N))/( JA(N)*CMPLX(DJ(N) 
C  ,-DY(N))  -  ZR*DJA(N)*CMPLX(J(N),-Y(N))) 

WRITEC*, 1000)  N,GAMMA(N.1),GAMMACN,2) 

10  CONTINUE 

C 

C  CALCULATING  SIGMAs 

C 

DO  30,  I  =  180,  -180,  -ARES 
DO  40,  II  =  1,  2 

SIGMA(II)  =  (0. 0,0.  0) 

DO  20,  N  =  1,  MODE 

SIGMA(II)  =  SIGMA(II)+2.  0*GAMMA(N,II)*COS(N*PI*I/180.0) 
20  CONTINUE 

SIGMA(II)  =  SIGMA(II)  +  GAMMA(0,II) 

SIGMAN(II)  =  ((4.0/K0)*(CABS(SIGMA(II)))**2) 

40  CONTINUE 

WRITEC 3, 120)  I,  SIGMANCl),  SIGMANC2) 

30  CONTINUE 

C 

C 

110  FORMATClX.’Er  =  ' ,F6. 4, 1X,F6. 4,/, 

C'  Mu  =  ',F6.4,1X,F6.4,/, 

C  RADIUS  CMETERs)  =  ',F8. 5,/, 

C  MAX  MODE  =  ' ,13,/, 

C  KO  =  ',F8.5,/, 

C  KrRa  =  ' ,F8. 5,1X,F8. 5) 

120  F0RMATC1X,I4,2C3X,E14.4)) 

1000  F0RMATC1X,I3,1X,2CE12.4,1X,E12.4,4X)) 

C 

C 

STOP 

END 
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APPENDIX  G.  SOFTWARE  SOURCES 

1.  DISSPLA 

Integrated  Software  Systems  Corporation 
10505  Sorrento  Valley  Road 
San  Diego,  CA  92121 

2.  CURVE-DIGITIZER 
West  Coast  Consultants 

4202  Genesee  Avenue,  Suite  309 
San  Diego,  CA  92117 

3.  Microsoft  FORTRAN 
16011  NE  36th  Wav 
BOX  97017 
Redmond,  WA  98073 

4.  Microwav  \DP  FORTRAN 
POB  79  ' 

Kingston,  MA  02364 

5.  Prof  M.A.  Morgan,  Code  62Mw 
Naval  Postgraduate  School 
Monterey,  CA  93943 

6.  LT  T.B.  Welch  III 
C  O  T.B.  Welch  Jr. 

1318  Walthour  Road 
Savannah,  GA  31410 
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